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SUMMARY OF WORK 


The application of large-scale simulation to problems in shape optimization is very difficult 
for (at least) three reasons. First, the objectives can not always be stated a priori, and thus 
there is reluctance to expend resources. Second, simulations are still very expensive, at 
least for most complex physical problems, and thus optimization — which calls for 
repeated appeal to the forward problem — is often prohibitively costly. Third, many 
interesting problems in shape optimization also require variations in topology; the latter 
introduces singularities that frustrate many classical optimization procedures. 

In this grant, we have proposed an approach to shape optimization that addresses, or at 
least mitigates, these difficulties. First, as regards how the problem is posed, we have 
applied concepts from multi-criterion optimization theory, in particular, Pareto theory. The 
Pareto formulation permits the designer to identify (monotonic) preferences in certain 
performance metrics (e.g., lift and drag), and then obtain a trade-off curve between these 
variables; all points on this trade-off curve are optimal in the sense that there is no design 
point at which both metrics can be improved (e.g., both lift increased and drag decreased). 
This trade-off curve can then be used to help the designer identify the right balance between 
different preferences, and determine the optimal operating point once these preferences are 
established. 

Pareto optimization has been used only rarely (if ever) for complex problems in fluid 
dynamics because the determination of the trade-off curve constitutes a difficult min-max 
problem that requires many appeals to the (expensive) simulation. Our approach to 
alleviating this computational bottleneck is to replace the simulation with a simulation 
surrogate. In particular, we first construct a model, or surrogate, of the input-output 
behavior of the simulation — based on some small set of input-output pairs — and then use 
this surrogate, not the original simulation, in the resulting (Pareto) optimization. The 
surrogate, unlike the original simulation, is very inexpensive to evaluate, and thus 
extensive optimization can be performed; of course, the surrogate may also be significantly 
less accurate than the originating simulation. 

To address the accuracy question, we have developed an extensive statistical validation 
procedure and associated theory for understanding how well a surrogate is performing, and 
whether any particular surrogate-predicted result can be trusted. The latter is, in fact, 
facilitated by the Pareto framework, since the region of input space of interest is narrowed 
to the pre-image of the trade-off curve. In particular, if the number of performance metrics 
of interest (e.g., lift, drag) is small, the Pareto pre-image will typically be a low- 
dimensional manifold, even if there are many design variables (inputs). 

Finally, in order to treat the topology issue, we have implemented a level-set geometry 
description within the surrogate context. In this approach, shapes are described as level 
sets of parametrized functions such that (discontinuous) changes in geometry' can be 
described by continuous changes in parameters which, in turn, correspond to continuous 
changes in the performance metrics. The latter thus permits the use of gradient 
information, critical to any efficient search algorithm. 

The ingredients described above have been demonstrated in the simulation-based 
optimization of a compact (laminar-flow) heat exchanger. The performance metrics in the 
Pareto analysis are the pumping power and heat removal; the simulation is a full unsteady 
Navier-Stoices calculation, on the basis of which a simple (radial-basis function) surrogate 
input-output model is constructed; level sets are used to describe the geometry of the eddy 
promoters through which the flow is excited and the thermal-hydraulic performance 
improved. 




The results of this study indicate that our approach can prove quite effective in practice, 
permitting us to address problems not amenable to other techniques. However simple 
("connect-the-dot") surrogate methods are fundamentally limited to rather low-dimensional 
inputs spaces; although the Pareto approach permits us to validate a surrogate over a low- 
dimensional manifold, we must first construct the surrogate over the entire input space. As 
is well-known, approximation in many dimensions is plagued by the "curse of 
dimensionality" — there are simply too few points in each coordinate to represent the 
function unless a prohibitively expensive (exponential) number of input-output pairs is 
used. 

One promising approach to improving the dimensionality scaling of surrogates is to exploit 
state-space based models — models that contain not only simple input-output information, 
but also information originating in the underlying mathematical system that describes the 
phenomenon. In a follow-on grant we will be studying this possibility, using for our 
"surrogate" a low-dimensional numerical approximation, and for our validation procedure a 
new form of a posteriori error estimation theory. The latter also replaces the statistical 
results of our earlier work with purely deterministic bounds (albeit at some loss in 
generality). 

A detailed description of the results of the current grant are included in Appendix A in the 
form of a recently completed Ph.D. thesis by John Otto. 
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Abstract 

A nonparametric- validated, surrogate approach to optimization has been applied to the 
computational optimization of eddy-promoter heat exchangers and to the experimental op- 
timization of a multielement airfoil. In addition to the baseline surrogate framework, a 
surrogate-Pareto framework has been applied to the two-criteria, eddy-promoter design 
problem. The Pareto analysis improves the predictability of the surrogate results, preserves 
generality, and provides a means to rapidly determine design trade-offs. Significant contri- 
butions have been made in the geometric description used for the eddy-promoter inclusions 
as well as to the surrogate framework itself. 

A level-set based, geometric description has been developed to define the shape of 
the eddy-promoter inclusions. The level-set technique allows for topology changes (from 
single body, eddy-promoter configurations to two-body configurations) without requiring 
any additional logic. The continuity of the output responses for input variations that 
cross the boundary between topologies has been demonstrated. Input-output continuity 
is required for the straightforward application of surrogate techniques in which simplified, 
interpolative models are fitted through a construction set of data. 

The surrogate framework developed previously has been extended in a number of ways. 
First, the formulation for a general, two-output, two-performance metric problem is pre- 
sented. Surrogates are constructed and validated for the outputs. The performance metrics 
can be functions of both outputs, as well as explicitly of the inputs, and serve to char- 
acterize the design preferences. By segregating the outputs and the performance metrics, 
an additional level of flexibility is provided to the designer. The validated outputs can be 
used in future design studies and the error estimates provided by the output validation 
step still apply, and require no additional appeals to the expensive analysis. Second, a 
candidate-based a posteriori error analysis capability has been developed which provides 
probabilistic error estimates on the true performance for a design randomly selected near 
the surrogate-predicted optimal design. 


Thesis Supervisor: Anthony T. Patera 
Title: Professor of Mechanical Engineering 
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Chapter 1 


Introduction 


1.1 Motivation 

The unprecedented advances in both computational and simulation capabilities have spurred 
a shift in focus from the development, improvement, and application of single-analysis tech- 
niques. to the incorporation of high-fidelity simulations into formal design optimization 
frameworks. Although formal optimization techniques have been widely used in structural 
design [34. 75], they have only recently been used regularly in fluids applications. Com- 
putational fluid dynamics simulations are typically very CPL -intensive, and when directly 
inserted into optimization routines as a function-call, the evolution of the design to a local 
optimal will advance extremely slowly. The lack of interactivity that characterizes such an 
approaches may be acceptable in the very late stages of the design, but is unacceptable in 
the early stages of design in which widely varying configurations are often examined and 
trade studies are typically conducted. 

Preliminary design is, intrinsically, an interactive process. It requires that a large num- 
ber of assessments be made rapidly as design parameters are changed and updated. Fur- 
thermore, as more is learned about the design in terms of its performance and reaction to 
changes in input quantities, goals are likely to evolve and change as well. Methods in which 
each design cycle requires hours, or even days, are clearly not useful during the evolutionary- 
stages of preliminary design in which flexibility and interactivity are most valuable. 

An approach that has been widely pursued in an attempt to shift as much information as 
possible from the late stages of the design process, where radical changes are much too costly 
to be practical, to the early stages, where many design scenarios and configurations are 
considered, is through analysis-approximation schemes. Fundamental to such approaches 
is the notion that the extremely high simulation cost can be front-loaded, and precisely- 
budgeted to the construction of simulation surrogates: inexpensive input-output functions 
that duplicate the responses of the high-fidelity simulation. The simulation surrogates 
can then be coupled with the formal optimization techniques, the result of which is a 
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(potentially) highly-accurate, interactive design framework. 

The primary drawback to surrogate-based optimization approaches is the uncertainty 
that is introduced, both in terms of how the performance predicted by the high-fidelity 
simulation differs from that of the surrogate near particular designs, predictability , and 
also as to how a surrogate-predicted optimal design point compares to the simulation- 
based optimal design point, optimality [7]. Beyond post-evaluation of surrogate-predicted 
design points with the simulation, both predictability and optimality are often times ignored 
in surrogate-based optimization methods reported elsewhere. The validation and error 
analysis steps of the nonparametric-validated surrogate framework described, and applied, 
in Chapters 4 — 5 seeks to address these issues. 

A second issue in design optimization is how to best parameterize the surfaces of the 
object to be optimized. A great deal of effort has focused on shape parameterizations. In the 
work presented here, the underlying desire to examine shape-families that include multiple 
topologies (single and two-body configurations), and to treat the entire family in a unified 
approach, has greatly restricted the types of methods that are appropriate. The problem 
central to this thesis is of an eddy-promoter heat exchanger. The desire to extend the 
design problems beyond the single cylinder designs examined in previous work [41, 77. 79]. 
and the potential beneficial interaction that may result from two-body configurations, has 
motivated the requirement that the shape-family encompass multiple topologies. 


1.2 Design Optimization 

Formal design optimization procedures have received increased attention in the fluid dy- 
namics community recently. Both the increased computer capabilities and the improve- 
ments in computational fluid dynamics codes in terms of robustness and capabilities have 
contributed to the gradual shift to formalized design optimization. In general, the optimiza- 
tion approaches that have been pursued can be broken into two broad categories: on-line , 
direct insertion approaches in which the CFD simulation serves as a function call inside of 
the optimization routine, and off-line , surrogate approaches in which the CFD simulation 
is used to construct simplified input-output models that are then used in the formal opti- 
mization routines. Each approach has associated advantages and drawbacks that make the 
choice between which is best problem dependent. 

The primary advantage to on-line optimization strategies is that they can handle prob- 
lems with a large number of input (design) variables very efficiently. Optimization methods 
are typically gradient based, and restricted to local searches, which reduces the effective 
dimension of the problem to that of the low-dimensional manifold traversed during the 
optimization process. Issues as to whether a local or global optimal point is achieved at 
the end of the process remain, but at the least, the local searches can be accomplished with 
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a moderate number of appeals to the analysis even for problems with very large (in terms 
of independent inputs) problems. An additional advantage, in the context of the compari- 
son between on-line and off-line strategies, is that upon completion of the on-line design 
optimization, a true optimal point is obtained. 

The biggest drawback to on-line optimization approaches is the inordinate amount of 
CPU time that they usually require. The underlying analysis codes are often very expen- 
sive, requiring on the order of hours of CPU time for a single analysis case. Not only is the 
computational time inordinately large, it is very difficult to know a priori how many sim- 
ulation evaluations will be required by the optimizer, making accurate resource budgeting 
impossible. Scenarios in which computational resources are exhausted before a meaning- 
ful design is obtained can be easily envisioned. Furthermore, as more is learned about 
the design during its evolution, goals and specifications often need to be modified as well. 
Changes in design parameters require an entirely new set of simulation analyses to perform 
the requisite design optimization update. The large computational cost of the on-line tech- 
niques greatly reduces the flexibility of the design process and makes interactivity out of 
the question. Secondarily, integration of the analysis code with the optimization routine 
can be difficult, especially for multidisciplinary optimization applications, irregularities in 
the analysis response can cause convergence problems, and finally, the incorporation of data 
from outside sources is difficult. 

Off-line, surrogate optimization approaches offer a number of advantages that directly 
address the disadvantages of on-line optimization. First, by design, the input-output sur- 
rogates are trivial to evaluate, ensuring that the optimization process can be completed in a 
reasonable amount of time and global optimal points can be obtained. The computational 
burden is placed on the construction of the input-output surrogates which occurs prior to 
any optimization. Because surrogate construction is performed off-line, and because the 
cost of a single analysis can be accurately estimated, available resources can be precisely 
budgeted. Design goals can be changed, and new optimal points obtained, throughout the 
process, and interactivity is preserved. Integration is greatly simplified as all information is 
contained in a simplified input-output function, smoothness can be enforced on the surro- 
gates, and data from outside sources is readily incorporated into the surrogate construction 
process. 

As already stated, the primary drawback to approximate optimization strategies is a 
consequence of the surrogate-for-simulation error. The construction of surrogates in high 
dimension input spaces is difficult, and as will be shown in Chapter 4. design predictability 
deteriorates rapidly with increased input dimension. The validation and a posteriori error 
analysis stages of the surrogate framework seek to address the errors that are incurred when 
optimizing with respect to a simplified model as opposed to the high-fidelity analysis. 
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1.2.1 Approximate Optimization Techniques 

Approximate optimization techniques have been pursued in wide range of optimization con- 
texts [3]. Central to the notion of approximate optimization is that a (typically) expensive 
analysis represents the truth. Direct insertion of the analysis code into the optimization 
routine is prohibited by the high cost associated with each appeal to the analysis and a sim- 
plified model (a surrogate) is typically used in place of the analysis. The only requirement is 
that the surrogate approximate the relevant input-output relationships of the analysis. An 
obvious approach is to sample the analysis at a set of input points, obtain the corresponding 
outputs, and form a response surface [10]. The response surface serves as a simulation sur- 
rogate in the optimization [21, 67]. The surrogates can also be lower fidelity, less expensive 
simulations. Variable-complexity, response surface strategies [32] incorporate a hierarchy 
of approximations into the optimization process. 

The distinction that is made between the high-fidelity analysis and the simplified ap- 
proximation can be made more general. For example, an experiment may serve as the truth 
and a simplified input-output model would be the surrogate. An application of this type 
has been examined in Chapter 6 and is also reported in [54]. 

1.2.2 Nonparametric-Validated Surrogate-Based Optimization 

The key feature that is lacking in many approximate optimization strategies is a rigorous, 
well-integrated means for error estimation. Post evaluation of the surrogate-predicted 
designs is widely used and provides the requisite information to determine how a design 
actually performs. However, many times even a single simulation evaluation will destroy 
the interactivity of the design process. In situations in which several designs are selected, 
the cost of appealing to the analysis for each design become prohibitive. Also, in many 
situations the analysis may not be readily available during the design phase (certainly the 
case for the experimental optimization problem in Chapter 6). Finally, while post evaluation 
does inform the designer as to how the selected point performs, it provides no guidance as 
to what should be done if design goals are not met and does not give any information as 
to how nearby points perform. Post evaluation does not address how surrogate-predicted 
optimal designs relate to true optimal designs. 

The nonparametric-validated surrogate framework seeks to address the surrogate-for- 
simulation error. The baseline framework [77, 79] can be broken into three steps: (1) 
Off-line surrogate construction and validation. (2) Surrogate-based design optimization. 
(3) A posteriori error analysis. The framework can accept and assess any combination of 
surrogate and truth simulation (or experiment), and is valid regardless of the approximation 
quality of the surrogate. The results are rigorous, based only on verifiable assumptions, and 
require no additional parameters which typically must be estimated. Additionally, the error 
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estimates apply to the specific design (or region near the design) and precisely quantify the 
response. By contrast, root-mean-square errors used in the context of response surfaces 
[10] provide an assessment of the overall surrogate approximation accuracy but can not be 
easily applied to analysis of surrogate-predicted optimal designs. 

he baseline surrogate framework was first developed by [77. 79] for non-noisy simula- 
tions. It has been extended to the analysis of simulations with noisy outputs [78, 77] in 
which the output error is symmetric and unbounded. The original, baseline framework 
has been extended to a more general, scaled-error formulation for multiple outputs [57] 
and more sophisticated elemental, sequential, and adaptive techniques have been developed 
as well [55. 56]. Finally, surrogate methods have been coupled with Pareto analysis for 
problems with multiple, competing performance goals [36. 37]. 

The surrogate methods have been applied to a wide range of problems. The eddy- 
promoter heat exchanger that is examined in Chapters 2 — 5 has also served as an illustrative 
problem in [36, 77, 79]. The baseline surrogate framework has been used to optimize the 
profile of axisymmetric bodies of revolution in Stokes flow [55. 56. 57]. The noisy output 
formulation has been applied to conduction in random media [78. 77] and the elemental and 
sequential/adaptive techniques have been applied to the conduction problem of a composite 
[55] and design of trapezoidal ducts [56]. 

In this thesis, the baseline surrogate framework with error scaling, given in [57] is first 
applied to the eddy-promoter heat exchanger design problem that is characterized by two 
performance metrics. The performance metrics are functions of outputs from the simulation 
as well as explicitly of one of the design inputs. The surrogate framework has been extended 
to assess a general two output, two performance metric problem in which the simulation 
outputs are validated. By validating the outputs instead of the performance metrics, an 
additional level of flexibility is afforded the designer in that the validated out put can be used 
in later design studies with modified performance metrics. A distinction is made between 
the inputs which are modeled to construct, the output surrogates, and those for which the 
performance metric response is known analytically. By separating the design inputs in 
this way, the greatest advantage can be taken of what is known analytically. The surrogate 
theory has also been extended to assess the error on a design point randomly drawn near the 
surrogate predicted optimal design. The original surrogate framework made an assessment 
of the surrogate error incurred for regions of points near the surrogate-predicted optimal 
point. 

Second, because the eddy-promoter problem is characterized by two. competing perfor- 
mance metrics, and because the output surrogates and hence the performance metrics are 
very inexpensive to evaluate, a Pareto formulation of the problem is pursued. An extremely 
valuable synergy between the Pareto analysis and the surrogate methods has been identified 
[36] in previous work. In high input dimensions, the predictability results of the surrogate 
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framework become very poor. To alleviate the loss in predictability, a Pareto analysis is 
pursued which reduces the effective dimension of the problem from the number of inputs 
to one less than the number of performance metrics. For the two performance metrics 
eddy-promoter problem, the Pareto analysis reduces the optimization problem to that of a 
single parameter and sacrifices very little generality in the process. A fine-grained Pareto 
analysis of the eddy-promoter problem using the simulation is totally out of the question 
and a surrogate strategy is necessary to make the problem tractable. 

Finally, in Chapter 6, the baseline surrogate framework is applied to the experimental 
optimization of the three-element airfoil. The experiment is assumed to be deterministic, 
and the non-noisy formulation is used. 

1.3 Shape Optimization 

The eddy-promoter, heat exchanger that serves as the central example for the majority 
of this thesis falls into the wide class of shape optimization problems [63, 65]. Shape 
optimization problems are present in most engineering disciplines. The general structural 
design problem of determining the optimal truss is one such problem. For the truss problem, 
the thickness, length, positions, and possibly the existence of the various components must 
be determined and serve as inputs. Numerous shape problems can be found in the field 
aerodynamics. The optimization of a wing (both planform and the cross-sectional airfoil 
shape) to achieve desired lift, drag, and stall characteristics is frequently considered and 
of critical importance for aircraft design. The design of a wind-tunnel nozzle contour to 
achieve a prescribed flow profile in the test section is another example of shape optimization. 

1.3.1 Topology 

An important consideration in some shape optimization problems is that of topology. In 
truss problems [6], not only do elements need to be sized in terms of thicknesses and lengths, 
but connectivities and existence must also be considered. This extends truss optimization 
from simply that of sizing to include topology of the truss components as well. Several 
approaches to topology optimization have been demonstrated on structural optimization 
problems. The first approach is the ground structure method in which a given set of 
connections is allowed between a fixed set of nodes, and the problem is posed as a continuous 
optimization problem in which sizes can go to zero. Zero sized truss components imply no 
connection between the two nodes [4]. This approach is restrictive in that the connectivity 
(or the family of allowable connectivities) is assumed at the outset. 

The second approach to topology optimization, homogenization methods [4, 5, 16, 72], 
seeks to treat the problem in a more unified, pseudo-continuous sense. In a homogenization 
formulation, the distribution of a material property (e.g. modulus or density) serves as the 
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design variable, and the design is optimized for a given objective. Regions of the structure 
in which the value of the material property falls below a given cutoff are assumed to be 
voids in the structure. Because an underlying mesh is required, the resulting truss must 
often times "cleaned up” to smooth boundaries. The post-optimization smoothing that is 
required calls into question the actual optimality of the resulting structure. 

A third strategy used to find the optimal topology of a structure, zero-one methods [1. 
29], discretizes the solution domain into a a mesh of points and each point is systematically 
turned "on” and “off” until an optimum is achieved. The problem is thus reduced to the 
combinatorial optimization problem of finding the best set of material points that satisfy 
the constraints and minimizes the objective. 

While such treatments of topology optimization have been effective in many instances for 
solid mechanics applications, they are not directly applicable to fluid problems and are still 
more unsuitable for surrogate-optimization approaches. Topology is an important consid- 
eration. however, for some classes of fluid problems and should therefore not be neglected. 
The particular application explored in this thesis, the eddy-promoter heat exchanger, is 
one example of a problem that can benefit greatly from a unified treatment of multiple 
topologies. In particular, situations in which two bodies may favorably interact, and which 
may be more efficient than single body configurations, can be easily envisioned. 

None of the strategies pursued for the structural topology optimization are directly rel- 
evant or applicable to the eddy-promoter problem. Homogenization techniques do not have 
an analogous extension to fluid problems that is physically reasonable. The discrete zero one 
methods could be applied but suffer from several drawbacks. First, for the very expensive 
analysis that is used for the eddy-promoter problem, the combinatorial optimization prob- 
lem would be prohibitively expensive. A sufficiently refined discretization of the solution 
domain would make the problem still more intractable. Second, the zero-one techniques 
are not at all applicable to the surrogate techniques (e.g. a smooth input-output function 
can not be fit through a 20 x 20 mesh of discrete zero-one inputs) that are required to make 
the problem computationally possible. Third, and finally, the resulting optimal shapes are 
non-smooth and are likely not optimal in the sense of viscous dissipation. Post-processing 
of the optimal shapes to smooth the surfaces could be used, but relationship between the 
smoothed bodies and the true optimizers is unknown. The incompatibilities between the 
structural topology methods described above and the eddy-promoter design problem have 
motivated the development of a level-set based geometry description. 

1.3.2 Level Set Geometry Description 

The approach that has been developed to describe the geometry sacrifices the generality of 
the homogenization and zero-one techniques, but is more relevant to the eddy-promoter 
class of fluid problems. It is a level-set based technique and has several advantages. First, 
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it provides an unified (albeit restricted) approach to topology optimization in that no addi- 
tional logic is needed to distinguish between single or multiple body configurations. Second, 
it is directly applicable to surrogate approaches as it is efficient in terms of the number of 
inputs that are required and, most importantly, the input-output function is continuous 
across a topology change. Third, the family of shapes that are described by the level-set 
method are a superset of the cylindrical family studied elsewhere [40, 77, 79] making the 
results directly relevant to previous studies. Finally, it can be extended to more complex 
geometries and topologies in a rational way. The geometry is defined using superposi- 
tion principles and a core, generating function. By systematically adding core functions to 
the overall shape function, richer families of bodies can be defined, and a wider class of 
topologies can be considered. For the problem examined here, two distinct topologies are 
considered: single-body and two-body configurations. 

Level-set methods have been used previously to describe the domain in shape opti- 
mization problems [12, 80]. More recently, level-set methods have been widely applied to 
problems that involve propagating surfaces or interfaces [53]. This includes problems such 
as flame propagation and material boundary evolution. Level set methods have also been 
applied to character recognition (the determination of true character boundaries) [52] and 
image enhancement in the presence of noise. One of the principle strengths of level-set 
methods is their ability to handle difficult geometries including sharp corners, cusps, and 
topological changes [13]. The specific method used to describe the geometries is given in 
Section 2.5. 

1.4 Overview of the Thesis 

In the first chapter of this thesis, the eddy-promoter problem is described. The configu- 
rations, governing equations, shapes descriptions and summary of the fundamental design 
inputs, outputs, and performance metrics are included. In the second chapter, the numer- 
ical solution techniques used to solve the governing equations are detailed. The baseline 
surrogate framework is applied to the eddy-promoter problem in the third chapter and the 
Pareto surrogate approach is applied to the eddy-promoter problem in the fourth chapter. 
At the end of the fourth chapter, key features of the baseline and Pareto approaches are 
compared and contrasted. In the fourth chapter, the baseline surrogate approach is applied 
to the experimental optimization of a three element airfoil. Finally, some concluding re- 
marks are given in the final chapter followed by appendices that contain surrogate proofs, 
details of the numerical implementation, and the surrogate input-output functions. 
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Chapter 2 


Eddy— Promoter Heat Exchanger 


A natural byproduct of almost all mechanical processes is unwanted heat. The amount of 
work in the field of heat transfer problems in general has reached such a critical mass that 
reviews of heat transfer reviews are available [19]. Additionally, heat exchangers are present, 
in almost all mechanical systems and have received a great deal of attention as well [25. 76]. 
The requirement that the heat be efficiently dissipated to preserve mechanical integrity 
without incurring undue performance penalties has helped to support a very active field of 
research focused on heat transfer enhancement. 

One approach to the removal of heat from a surface is through a laminar flow heat 
exchanger. In such an exchanger, heat is transfered from a generation source into a fluid 
medium and carried downstream. A number of ways exist to enhance the heat transfer away 
from the wall of the channel and into the fluid medium. These include active mechanisms 
such as oscillatory pumping of the fluid or passive enhancement in which wall contouring or 
flow obstacles are introduced into the fluid channel to increase mixing. For all of the meth- 
ods. the intent is to increase the convective heat transfer by increasing the mass transport 
away from the wall. 

The heat exchanger geometry introduced in this chapter, and used as an example in 
Chapters 4 and 5. uses a passive enhancement strategy in which flow obstacles are intro- 
duced to lower the critical Reynolds number. The exchanger, as defined here, could be a 
subsystem in a much larger industrial heat exchanger device or could be a very close approx- 
imation to a design that might be used in a space-limited application such as the cooling 
of electronic devices. The enhancement geometries studied here take as a departure point 
work performed earlier in which cylindrical inclusions (eddy-promoters) were studied [41]. 
The cylindrical eddy-promoters greatly reduce the Reynolds number at which flow instabil- 
ities develop. The onset of instabilities have been observed at Reynolds numbers as low as 
125 when the eddy-promoters are present [40]. The presence of the eddy-promoter inclu- 
sions greatly increases the temperature convection in the channel and therefore, improves 
the heat transfer. The family of eddy-promoter geometries studied here include cylindrical 
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Figure 2-1: General geometry for the heat-exchanger configuration. 


inclusions, oblong inclusions at various orientations, and two body configurations. A uni- 
fied geometry description has been developed to describe the full family of geometries that 
preserves continuity in the input-output functions. The geometry description is based on 
level-set methods and is described in Section 2.5. 

In this Chapter, the general configuration and design goals for the eddy-promoter heat 
exchanger problem are presented. A description of the heat exchanger configuration and 
the simplifications made for later computations are presented first. Second, the governing 
equations and the scalings that are made to reduce the problem to nondimensional form 
are described. Third, the engineering goals for the problem are discussed and the nondi- 
inensional formulation for the design problem are given. Fourth, the Pareto formulation is 
presented followed by the level-set based geometry description. Finally, the reduced for- 
mulation for the eddy-promoter design problem is presented in terms of the design inputs, 
design spaces, outputs, and performance metrics. 

2.1 Physical Problem 

The general configuration of the heat exchanger is given in Figure 2-1. The flow is from left 
to right and is driven by an applied pressure gradient — The channel length is V and 
the channel half height is H 9 . The temperature of the fluid at the inflow is T/ n , the wall 
temperature at the outlet is T' w and the channel average velocity V f is 

1 f 2H ‘ 

v =wL u ' iy ’ (21) 

where u' is the x'-component of the vector velocity u' = (u',r/). There is is a uniform heat 
flux q n through the lower wall and the upper wall is assumed to be insulated. The objective 
is to transfer as much heat as possible through the lower wall, into the fluid stream, and 
out of the heat exchanger channel. 

A number of methods exist to effect better transfer of heat into the fluid medium. These 
include unsteady pumping, grooved walls, and the inclusion of periodically-spaced, cylin- 
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Figure 2-2: General geometry eddy-promoter, heat -exchanger configuration (top) showing 
one periodicity cell (bottom). 


drical obstructions referred to as eddy-promoters. In this work, the focus will be exclusively 
on the eddy-promoter configurations, a segment of which is shown at the top of Figure 2-2. 
The periodicity spacing (the length between inclusions) for the eddy-promoters is l' . The 
primary effect that the eddy-promoters have on the flow is to trip instabilities at greatly 
reduced Reynolds numbers than observed for plane channel flows. Two-dimensional, lin- 
ear instabilities are observed in channel flows at Reynolds numbers based on the channel 
half height H' of approximately 5500 [17, 60]. For properly placed and sized cylindrical 
eddy-promoters, the onset of two-dimensional flow instabilities have been observed numer- 
ically and experimentally at Reynolds numbers as low as 125 [40. 41]. The unsteady flow 
increases the convective transfer of heat away from the wall and increases the efficiency of 
the exchanger. 

As already stated, the heat transfer can be increased by the inclusion of obstacles in the 
flow to induce instabilities at lower Reynolds numbers. While the periodically spaced flow 
obstacles do increase the heat transfer, the increase comes at the cost of a higher pressure 
drop and greater viscous dissipation. The result of the higher pressure drop is that a more 
powerful fluid pump must be employed, increasing the cost and complexity of the exchanger. 
Additionally, excessive viscous dissipation, pressure forces, and shear stresses produce larger 
structural loads that must be balanced by stronger materials to avoid fatigue failures. 

The effectiveness of the exchanger can be quantified by two characteristics; how much 
heat it can transfer away from the hot wall and how much effort is required to accomplish the 
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transfer. The two quantities that characterize the exchanger are, in most instances, at odds 
with each other. Clearly, the heat exchanger design is an example of a thermal-hydraulic 
trade-off problem. In words, the goal of the heat-exchanger problem is to maximize the 
heat transfered into the fluid medium and carried downstream while holding the required 
pumping power as low as possible. These competing goals form the trade-off problem. 


2.2 Governing Equations 

The total exchanger length L' is assumed to be such that the flow is fully developed in 
x'. Additionally, the flow is assumed to be independent of the spanwise coordinate and is 
modeled as a two-dimensional problem.. Furthermore, for a periodicity spacing, l' , between 
eddy-promoter inclusions that is also sufficiently long, the fluid solution will be also be l'- 
periodic. With these two assumptions, the full exchanger length can be modeled with a 
single periodicity cell that enforces periodic boundary conditions at each end. 

The flow solution for the full channel shown in Figure 2-2 is modeled by a single peri- 
odicity cell as depicted in the lower half of Figure 2-2. The fluid flow solution is governed 
by the incompressible Navier-Stokes equations and it is assumed that all fluid properties 
are held constant and are independent of the fluid temperature. The equations governing 
the fluid flow, in dimensional form are 

+ pu' - V'u' - pV'V + VV - 0, (2.2) 

V' u' = 0. (2.3) 

where p is the fluid density, p is the fluid viscosity, and n' is the pressure. No-slip boundary 
conditions are enforced on the velocity along solid walls 

u' = 0 on dVf p \ dVf ", dVg”, (2.4) 

and left-right periodicity is enforced on the flow velocity and derivatives 

u'(x' = 0, y, t ') = u'(x' = /', y', t'), (2.5) 

K'i x ' = 0,y',t’) = vl x ,{x' = l',y',t'), ( 2 . 6 ) 

n' y ,{x = 0,y',t') = u ' y ,(x' = l',y',t'). (2.7) 

The pressure tt'(x', t') has the form 

AP' 

t')=p'(x',t')-^x', (2.8) 
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where x' = (x 1 . y') is the vector of coordinate directions in two-dimensions and is the 

pressure gradient that drives the flow. The periodic part of the pressure, p'(x'.f'), satisfies 

p'{x = 0. y . t ') — p(x' — l. y' . t'). (2.9) 

The temperature is governed by the scalar, convection-diffusion equation 

+ V' • (u'T') = aV 2 T'. (2.10) 

at' 

where T' is the full temperature solution and a = is the thermal diffusivity. For this 
equation to hold, it is assumed that the thermal conductivity, k, is independent of the 
temperature, that diffusion can be neglected, and that there is no internal heat generation. 
The boundary conditions for the temperature are the prescribed uniform heat flux through 
the lower wall 

nV'T' ■ n = q" on dVf p ' . (2.11) 

and zero heat flux (insulating surfaces) through the eddy-promoter body surface and the 
upper wall 

kV'T' • n = 0 on <TD§ P '. <9X>P'. (2.12) 

The heat flux q" is assumed positive for heat transfer into the domain and n is the unit 
outward normal. 

The full temperature. T . is decomposed into two components 

T'(x'.t') = 0'(x',O+ 7 V. (2.13) 

The linear part of Equation 2.13 is the most general form for which there is a unique solution 
such that d'(x’.t') satisfies periodicity [30] 

0'{x' = 0 ,y',t') = 9'{x' = l,y',t'). (2.14) 

The derivation of the coefficient 7 ' is given below. 

The governing equation for 9' is obtained by substituting Equation 2.13 into the tem- 
perature equation 2.10 

^ + V' • (u'0') = aV 2 9' - Vu'e!, (2.15) 

at' 

where §1 = (1.0) is the unit vector in the x-direction. The boundary conditions on 9' are 
obtained similarly and are 

kV' 9' -n = q" on dvf 9 ' , 
kV9' ■ n = — «7 'n x on dV^' ■ dV ^ p> , 
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(2.16) 

(2.17) 



where n i is the component of the unit outward normal vector, n = (ni,n 2 ), in the re- 
direction. 

Because solution domain V EP ' is a single periodicity cell of the full heat exchanger shown 
at the top of Figure 2 - 2 , the coefficient 7 ' in Equation 2.13 must be set such that heat does 
not build up in the channel, on average, over time scales of the order of period of the 
steady-periodic state. To determine the constant coefficient 7 ', the temperature Equation 
2.10 is integrated over the solution domain V Epf 

f 0'dx'=-[ V'-[u'( 0 ' + 7 ' 'x)\dx'+[ aV 2 e'dx'. (2.18) 

dt f Jv EP ' Jv Ep/ Jv Ep/ 

The first integral on the right hand side of Equation 2.18 is expanded, and Green’s theorem 
is applied to obtain 


£f 

dt 9 Jv EP ' 


e f dx f = - f 

JdV EP ‘ 


u 'O ' « n dl f — 7 ' [ u f dx! 
Jv Ept 


L 


dV Epl 


aV'e'-ndx'. 


(2.19) 


The requirement that T' be bounded as t — > oc is equivalent to enforcing that the left 
hand side of Equation 2.19 be zero when averaged over a length of time on the order of the 
solution period. The first term on the right hand side of Equation 2.19 goes to zero by the 
velocity boundary condition given in Equation 2.4 and the second term on the right hand 
side can be rewritten in terms of the flow rate 



( 2 . 20 ) 


For the problem examined here, the flow rate is fixed and Q'{t) — Q' . The expression for 
j' is obtained by requiring that the remaining terms satisfy Equation 2.19 which can be 
rewritten as 

j'l'Q' — a f Vd'-ndl'. (2.21) 

J dv EP ' v 

The boundary condition contributions from the body and upper wall in Equation 2.17 are 
substituted into the right hand side of Equation 2.21. The contribution from the body is 
zero because of integration around the closed contour dV^ p and the contribution along the 
upper wall is zero because n\\ dv ep = 0. The boundary condition contribution from the 
lower wall in Equation 2.16 is substituted into Equation 2.21 to obtain 


= a f Q -dl , 

JdVf pf K 


( 2 . 22 ) 


which is solved for 7 ' to get 

7 = 2 ac Wv r ( 2 ‘ 23 ) 

where the relation Q' = H'V' has been substituted, and q' = I'q" is the integrated heat flux 
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into the domain V EP through the lower wall dT>f p . 

To this point, all quantities and equations have been given in dimensional form. To 
nondimensionalize the governing equations, and to simplify and generalize the solution 
process, the following scales are used. All lengths are scaled by the dimensional channel 
half-height H ' resulting in a nondimensional channel height of H = 2. The velocities are 
scaled by \V ' which has been selected so that the nondimensional centerline velocity for 
plane Poiseuille flow is unity. The dimensional time is scaled by the ratio of the length scale 
to the velocity scale and the temperature is scaled by 2 -^-. The scales are summarized 
below. 

Velocity scale: |V r ' 

Length scale: H ' 

Time scale: | yr 

tt 

Temperature scale: ^ — 


The relative periodicity lengt h /' and the full channel length V are scaled by H' and resulting 
nondimensional lengths are 


l = 


l' 

H' 


L = 


L' 

H' 


(2.24) 


respectively. 

The scalings listed above are applied to the governing equations 2.2. 2.3. 2.15 and to 
the associated boundary conditions 2.4 — 2.7. 2.14. 2.16. and 2.17. The resulting equations 
are simplified to obtain the full set of nondimensional governing equations and boundary 
conditions. The velocity and pressure are governed by the incompressible. Navier-Stokes 
and continuity equations 


?+u Vu- -t^V 2 u + Vp = f (<). 
at Re 

V • u = 0. 


(2.25) 

(2.26) 


with boundary conditions 


u = 0 on 


dV, 


EP 


ov 


EP 


0Vl P . 


(2.27) 


u(x = 0. y, t) = u(x = /. y , t). 
u x (i = O.y.f) = u x (i = l.y.t). 
u y (x = 0, y. t) = u y {x - l.y.t). 


The governing equation for the temperature reduces to 

1 


! + u.(V») 


RePr 


V~9 - ■yui, 


(2.28) 

(2.29) 

(2.30) 

(2.31) 
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with boundary conditions 


ve- 

n = 1 - 

7nj on dV pp , 

(2.32) 

V0-n 

= ~7«l 

on dVf, dV pp , 

(2.33) 

9(x 

= 0 ,y,t) 

= 9(x = 

(2.34) 

e x (x 

= 0 ,y,t) 

II 

H 

H 

II 

(2.35) 

0y(x 

= 0 ,y,t) 

Pi 

II 

II 

(2.36) 


The Reynolds and Prantl numbers are 

3 V'H' 


Re = 


2v 


Pr = 


a 


(2.37) 


and the Prantl number is held fixed at unity. Pr = 1, for all cases presented here. The 
nondimensional velocity, pressure, and temperature are 


u = 


2u' 

3V 77 


(H 


2 1 


and 9 — 


9 ' 

ijjnT' 


(2.38) 


respectively and the nondimensional coefficient. 7. becomes 

_ 3 

^ ARePr 


(2.39) 


The dimensional pressure n f has been expanded according to Equation 2.8. This yields 
a time-dependent forcing term on the right hand side of the ^-component of Equation 2.25 


AF' 



(2.40) 


where the full vector forcing term in Equation 2.25 is f (t) = (f(t).O). The forcing term is the 
nondimensional pressure drop required to drive the flow at the prescribed Reynolds number. 
To be consistent with the velocity scaling \ V , the nondimensional. time dependent, forcing 
term f(t) is set such that nondimensional flowrate given by 


=/ 

Jo 


u{x,y : t) dy, 


(2.41) 


is a constant Q = The numerical implementation that has been used to to achieve the 
prescribed flow rate is detailed in Section 3.2.1. 
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2.3 Engineering Goals 


The eddy-promoter heat exchanger design problem can be classified as a multicriteria design 
problem. The two criteria that characterize the design goals are the pumping power required 
to maintain the prescribed flowrate and the temperature rise along the lower wall. The 
two design criteria for this problem result in a classic trade-off problem in that design 
improvements, reduced pumping power and lower temperature rise, are competing goals 
and in many instances run counter to each other [66]. From the designer’s point of view, 
it is best to postpone how to precisely balance and prioritize the two criterion until as late 
in the design process as possible. This will increase the flexibility and allow the designer to 
learn as much as possible before making decisions. 

To accomplish the design goals, a number of independent variables can be adjusted. The 
independent variables are referred to as design (or input) variables. Of the full set of design 
variables, a subset consists of the vector of parameters Z' EP that specify the geometric 
configuration of the eddy-promoter inclusions shown shaded in the lower part of Figure 
2-2. The vector Z' EP includes the parameters that describe the profile of the body surface 
as well as its placement in the channel. Additional design variables are the height of the 
channel configuration H' . and the average fluid velocity V' defined in Equation 2.1. The 
total channel length L' is fixed as it is assumed to be set based on the specific application. 
The relative spacing between each eddy- promoter. can also serve as a design input but. 
for the remainder of the work presented here, it is also assumed fixed and is set such that 
1=4 'jj — 6.666. The value l = 6.666 matches that used in previous studies and has been 
shown to be a valid length for the existence of ^-periodic solutions [40]. The fluid properties 
p. v, k, and a are all assumed fixed. The full, dimensional design input vector is 

Z' = (Z ' ep .V',H'). (2.42) 

As stated in the first paragraph of this section, the two criteria that characterize the 
design are the dimensional pumping power required to achieve a prescribed average flow 
velocity V'' and the temperature rise along the lower wall. The pumping power is given as 

T' P (Z') = AP’V'H', (2.43) 


and is dependent on the input design vector Z. The temperature rise is broken into two 
components and is 


ST' 


= T' 

A W,0 


- T' 


(2.44) 


where T' w 0 is the temperature at the outlet of the heat exchanger, and T' m i = T' m (x = 0) is 
the mixed mean temperature at the inlet of the exchanger. The mixed mean temperature 



is defined as 

1 r 2H 

T' m {x) = J o < u'(x. t)T'(x, t) > dy. (2.45) 

The brackets < ■ > in Equation 2.45 indicate a temporal average over a sufficiently long 
period of time to obtain the mean of the time-varying, steady-periodic solutions. 

To nondimensionalize the problem, it is important to scale the optimization objectives 
by quantities that do not change as the design changes. This requirement simplifies the 
solution process but it also restricts the possible nondimensionalizations. The independent 
variables are scaled as well. The channel half-height H f is scaled by the total channel length 
L\ which is fixed. For convenience, the inverse of the scaled channel half-height 


V 



(2.46) 


serves as the design variable. The eddy-promoter configuration variables are scaled relative 
to the channel half-height 


Z EP = 


Z 'ep 

H' ' 


(2.47) 


The average velocity V' is scaled by jjj - which gives the Reynolds number. Re = . as 

a design variable. The full set of nondimensional design variables is then given as 


P full = (Z ep- Re.ru). 


(2.48) 


To scale the temperature objective, the temperature rise ST' is first expanded as 

q" l' 


ST' = 


pCpV'H' 


+ AT 


2.49) 


where AT ~ < T' m (x) — T^.(x) > and the over-bar indicates a spatial average in the x- 
direction over the length of the spacing between eddy-promoters, l'. The first term on the 
right hand side of Equation 2.49 corresponds to the mixed mean temperature rise along 
the exchanger length due to the uniform heat flux q" and the second term is the heat that 
has been transfered into the fluid medium. Scaling the temperature rise by gives the 
desired result for the objective 


0(Z ep- Re, L ) = 


ST' 


1 


+ 


1 


2^1 RePr t) L Nu{Z EP . Re) 

where the Nusselt number is given by 

Nu(Z E p, Re) = 


(2.50) 


1 


1 


(2.51) 


< T m (x) - T w {x) > AT 

and is the inverse of the nondimensional, temporally- and spatially-averaged difference 
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between the mixed-mean and wall temperatures. 

3 

The pumping power, T' P { Z'), is scaled by ^ 77 . The nondimensional pumping power is 

AP'V'H'L ' 2 

*(Z EP-Re,VL) = 3 • (2-52) 

pw 5 

Introducing the nondimensional friction factor 

AP' 

V’oIZep, Re) = — — -j, (2.53) 

and substituting the expression into Equation 2.52 gives 

*(Z EP. Re, tjl) = MZep- ReW L Re\ (2.54) 

where the nondimensional friction factor in Equation 2.53 is identically equal to the tempo- 
rally averaged, pressure forcing term on the right hand side of Equation 2.25. The pressure 
forcing term is given in Equation 2.40 and the output form is 

ipoC^EP- R e ) =< R e ) > ■ (2.55) 

To this point, no preference as to how to treat the two design criteria. 0(Z EP-Re.ru) 
and l P(Z ep- R e- Vi), has been stated. The precise application and engineering constraints 
(e.g. available pumps, limits on temperature rise, etc.) would need to be known to de- 
termine the best formulation for the design problem. As an example, if the maximum 
temperature rise along the lower wall of the exchanger must be limited (say due to oper- 
ating temperature limits on electrical components), the problem could be formulated as a 
constrained optimization problem 

min ^(Zep- Re.rji) (2.56) 

{Z E P ,Re,rj L } 

subject to 

0(Z E p.Re, VL ) <9, (2.57) 

where 9 is the upper-bound constraint on the temperature performance criterion. Con- 
versely, the problem could be formulated based on the maximum pumping power that can 
be supplied by the fluid pump. In that case, the two performance metrics in Equations 2.56 
and 2.57 would be reversed and the upper-bound constraint in Equation 2.57 would be set 
based on the particular fluid pump. Ideally, it is best to put off such decisions until as late 
in the design process as possible to preserve generality. The design strategy that is used 
for the eddy-promoter problem combines the formulation presented in Section 2.4 and the 
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surrogate techniques presented in Chapters 4 and 5. The result is an extremely flexible de- 
sign environment that addresses both constrained optimization problems mentioned above 
as well as a wide class of more general formulations. 

Regardless of the problem formulation that is ultimately selected, upon solution of the 
nondimensional optimization problem, the physical heat exchanger quantities must be ex- 
tracted from the nondimensional results. The steps are as follows: 1) Given the dimensional 
length of the heat-exchanger channel L\ H ' is determined from the optimization solution 
for t\l . 2) With H f known, the dimensional geometric configuration parameters Z f EP can 
be obtained. 3) The Reynolds number implies (for fixed fluid quantities) the optimal V 9 
given H' . 


2.4 Optimization Problem — Pareto Formulation 

The notion of Pareto-Optimality (PO) has its origins in the field of economics. The analysis 
was originally targeted towards the equilibrium conditions in an economy in which gains 
made by one consumer would be at the expense of at least one other consumer. The 
equilibrium point was first referred to by Edgeworth in 1881 [20] and the concept was 
refined and expanded by Pareto [59]. The application of Pareto analysis to engineering 
problems has only occurred recently, however. 

There are several advantages to pursuing a Pareto analysis of multicriteria optimization 
problems when feasible. First, it can be shown that a wide class of multicriteria optimization 
formulations have a solution that is also PO. This is important in that by pursuing a Pareto 
analysis, one loses very little generality. The second advantage to a Pareto approach lies 
in the engineering relevance of the results. For a two-criteria optimization problem, the 
resulting family of PO solutions map to a curve in the two-dimensional criterion-space that 
provides a designer with the capability to rapidly assess trade-offs and to determine the 
impact of changing design goals interactively. Third. PO analysis has a beneficial impact on 
the practicality of the error estimates provided by the surrogate framework. This is a direct 
result of the reduction in effective dimension that is realized by the Pareto analysis. The 
result is that for a problem with a high input dimension (M), the dimension over which the 
surrogates are validated is reduced to one less than the number of performance metrics that 
characterize the design preferences, (K — 1). The input dimension reduction significantly 
improves the predictability and is discussed in Chapter 4. 

In this section, the general formulation for the Pareto-optimization problem is presented 
first. Second, several classes of multiobjective problems that have solutions that are also 
Pareto-optimal are presented and discussed. Third, the methods used to determine the full 
set of Pareto-optimal solutions are described. Fourth, and finally, the generalization of the 
Pareto formulation to an arbitrary number of objectives is briefly presented. 
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2.4.1 Pareto-Optimal Solutions 


Central to the Pareto analysis of a problem is the notion that the designer can articulate 
monotonic preferences on each performance metric. For all of the discussion of this section, 
it will be assumed that there are two ( K = 2) performance metrics for which lower values 
correspond to better designs. Given a vector of design variables p € ft and the design 
(or input) space ft C 2R A/ , the two performance metrics are 4>i(p) = $i(<£i(p), <fe(p), P) : 
ft — ► IR and <f> 2 (p) = $l( < £l(p)t < fc(p);P) : ft -t 1R. It is arbitrarily assumed that the 
performance metrics are functions of the two outputs, </>i(p) : ft — > IR and <^> 2 (p) '• ft -K, 
as well as explicit functions of the input variables p. 

To define the Pareto-optimal sets, the notion of an achievable set 

A - {s G R 2 | 3 p € ft s.t. <M<Mp),<Mp):P) < •si-$2(<Mp)-<iMp)-P) < •‘ft}- (2.58) 

is first introduced. The achievable set A is the the set of all possible output pairs s = (si - s- 2 ) 
for which there is an input point p 6 ft that has correspondingly better performance metric 
values. With the achievable set defined as above, the set of Pareto-optimal designs is 
then the set of designs for which the performance metric values lie on the non-horizontal, 
non-vertical boundary of A not at infinity. This boundary, denoted dA. is depicted on 
the left in Figure 2-3 as the solid line and is termed the Pareto-optimal output manifold. 
The PO output manifold is also sometimes referred to as the efficient frontier [71] or the 
trade-off curve [66]. The inverse map of the PO output manifold back to the input space 
is determined from the performance metrics 

C A = {p € ft | 3 s € dA s.t. <Fi(tMp).<fe(p);P) = si.$ 2 (</>i(p)-</>>(p)-P) = «.’}• (2-59) 

and C A C ft is denoted the PO input manifold. The PO input manifold C A is shown as 
solid lines on the right side of Figure 2-3 and the outline of the design space ft is shown as 
a dashed line. Figure 2-3 depicts dA and C A as very well behaved. ( K - l)-dimensional 
manifolds. There is no justification for such behavior in general, especially as regards C A . 
but for relatively smooth, well defined performance metrics it is reasonable to expect some 
regularity. 

The design optimization interest is in the designs described by input points, p*. that lie 
on C A . p* 6 C A . These design points can rightfully be called the "best" or most efficient 
designs based on the performance metrics. This is illustrated by considering the input point 
p' with the corresponding output pair ($i (p'), $ 2 (p')) £ dA. For such a point, there always 
exists a point p* e C A that has output pair ($i (p*), <F>(p’)) S dA and that satisfies either 

•IMp*) < and $ 2 (p*) < 4> 2 (p') : (2.60) 
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Figure 2-3: A schematic the Pareto-optimal output (left) and input (right) manifolds. 


or 

$i(p*) < $i(p) and $- 2 (p*) < (2.61) 

In words. Equations 2.60 and 2.61 indicate that for a point p' that is not PO. there will 
always exist a point p* that performs at least as well in terms of the first performance 
metric and improves upon second performance metric <J> 2 (p 7 ) or (conversely) that 

performs at least as well as $2(p / ) but improves upon 5>i(p'). 

Finally, it is noted that the entire boundary of A does not correspond to PO designs. 
From the discussion of the previous paragraph, it is obvious that designs that correspond 
to horizontal or vertical segments of the boundary are not PO. The reason for this is that 
for such designs, there will always be another design input point that improves on one of 
the performance metrics without increasing the other. For all of the discussions of dA. and 
the corresponding C A , it is assumed that all points that are not PO have been eliminated. 
This has been suggested on the left side of Figure 2-3 by denoting dA as the solid line only. 


2.4.2 Relationship to Other Multi-Objective Problems 

The advantage of pursuing a Pareto analysis is that it is, in fact, a generalization of many 
multi-objective problems in that solutions to the given multi-objective optimization prob- 
lem, Pa /O’ are Pareto optimal: Pa/o £ C A . For example, the single performance metric 
minimization problems, 

p* = argmin$](p), (2.62) 

and 

P 2 = argmin$ 2 (p), (2.63) 
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both have solutions that lie on C A , p* € C A and € C A . In fact, the performance metric 
pairs that correspond to p] and P 2 lie at the extreme ends of .4. 

As a second example, if the objective is the weighted sum of the performance metrics, 

y?(p) = ai<£»i(p) + a2$2(p)j then 


p* = argmin<p(p) (2-64) 

pen 

is PO as well and can be found on the PO input manifold C A . p* 6 C A . 

In general, for any non-decreasing functions <£>(<I>i(p), $ 2 (p)r p) and 6 ($i(p)> $2(p): P) 
that are not explicitly functions of the input vector p. a solution to 

p* = argmin<pi($i(p),$2(p).p) ( 2 - 65 ) 

subject to 

6(<I>i (p), $2(p): P) < 0 (2.66) 

will lie on the PO input surface C A . This generality in the PO solutions allows the designer 
to eliminate uninteresting designs but ensures that sufficient flexibility remains to change 
the design goals as the design itself evolves. This flexibility is extremely important because 
many times, precise goals can not be defined at the outset of the design process. 

2.4.3 Determination of dA and C A 

To find C A . the multicriteria problem is reduced to a series of scalar problems parameterized 
by w G VV = [0. 1]. A min — max formulation for the scalar sub-optimization problems is 
used as it is guarantees that non-convex portions of C A are obtained [22, 14]. 

Given the scalarization parameter w G W, the scalar problem 

£(w) = arg min max^^i (p), (1 — tt;)4>2(p)), w € W, (2-67) 

pen 

is solved for a sufficiently refined W. The curve £(W) (with segments corresponding hor- 
izontal and vertical segments of the boundary of A removed) is the PO input surface C A . 
The corresponding PO output manifold is obtained from £(VV) and the performance metrics, 
given by 

OA = ($i( s c (W)), $ 2 (((W))). (2.68) 

It is immediately obvious that it can be an extremely time consuming and difficult pro- 
cess to find £(W) to an adequately fine resolution for even moderately complex (expensive) 
performance metrics. 4>i(p) and $ 2 (p)- Each solution to 2.67 is itself an optimization prob- 
lem. The prohibitive cost of obtaining £(W) is the primary reason why Pareto analysis is 
not typically pursued for complex engineering problems. 
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Note that in this discussion, it is assumed that the £(W) obtained from the scalarization 
process (2.67) has been “cleaned up” to remove duplicate points and horizontal and vertical 
segments of the boundary of A that, by definition, are not PO. 


2.4.4 Generalization to Arbitrary Output Dimension 

The notion of Pareto-Optimal solutions can be extended to an arbitrary number of per- 
formance metrics $i(p), • • • > $A'(p)- For the K performance metrics (for all of which it is 
assumed that lower is better) the PO output manifold is the {K - l)-dimensional manifold 
that represents the boundary of the achievable set Ak not at infinity. Here Ak is 

a k = {s e iR h | 3 p e ft s.t. <M0i(p),&(p).p) < si,... .$a-(<Mp),<£ 2 (p),p) < s K }. 

(2.69) 

In three dimensions, this can be visualized as the surface of the three-dimensional achievable 
set volume in the first quadrant that is “nearest" the origin. 

In higher dimensions, the surfaces can not be easily visualized but can be described 
with the min-max scalarization process similar to 2.67 [35]. Given the K — 1 scalarization 
parameters w = (uq, . . . , wk-i) € Wa _i, where 


the PO input manifold can be found by solving 


£(w) = arg min max ^^(p). . . . , u>a--i$a'-i(p)- ^1 - w,j 4>a'(p)J . w € Wa'-i, 

(2.71) 

where the scalarization vector space is Wa_i C [0. 1] 7 ' -1 . As before, the PO output 
manifold is given by 

dA K = (<M£(Wa-— i)). 4> 2 (^(Wa--i))). (2.72) 


and finite-sized portions of 8Ak for which at least one of the performance metrics is constant 
have been assumed removed. 

In practice, assuming very inexpensive performance metric functions, it may be possible 
to obtain and manage the data required to visualize and sample a three criteria ( K = 3) 
Pareto analysis. As noted above, for a three criteria optimization problem, the PO output 
manifold is a surface of a three-dimensional volume (the achievable set) and the PO input 
manifold can be visualized versus the two scalarization parameters as surfaces. With a 
sufficient amount of data, boundaries between disconnected subsets of C A can be accurately 
determined. For problems in which K > 3, the PO manifolds become very difficult to 
visualize (volumetric manifolds in four-space for K — 4) and dimensionalization issues make 
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it computationally difficult to precisely determine the boundaries between disconnected 
subsets of C A . With these considerations in mind, it can be concluded that the Pareto 
analysis is best suited for K = 2 problems and, possibly, could be useful for K = 3 problems. 
Fortunately, a large number of problems of interest can be meaningfully expressed as two- 
criteria optimization problems. 

2.5 Eddy-Promoter Shape Description 

The effectiveness of cylindrical eddy-promoters as heat transfer enhancement devices has 
been studied extensively [40. 77, 79]. For the eddy-promoter optimization problem studied 
here, the intent is to expand upon the the cylindrical family of designs that have been stud- 
ied previously. It is important that the definition of the design space include cylindrical 
geometries so that improvements to optimal cylindrical configurations can be identified. An- 
other desire is for the design space to include topologically changing geometries: specifically 
it should include single bodies as well as two-body configurations. 

The motivation behind examining two-body configurations is the possibility that the 
two bodies could favorably interact to improve the heat transfer beyond what is possible 
with a single body. The multiple-topology design space presents a particular challenge in 
that the definition of the shape parameterization must be consistent and continuous with 
respect to the outputs so that surrogates can be constructed. 

In this work, a shape parameterization has been developed that describes a family 
of geometries that is a superset of the cylindrical family studied previously. The shape 
parameterization presented encompasses a richer set of single-body geometries than simply 
circular cylinders, and includes two-bodv geometries as well. To accomplish this, a level-set 
based shape description has been developed. 

2.5.1 Level Set Methods - Shape Generation Example 

To describe the geometry of the eddy promoters, Z £p, a level-set [80, 64] based approach has 
been developed. Level-set methods have been applied to problems that involve propagating 
surfaces or interfaces [53]. This includes problems such as flame propagation and material 
boundary evolution. Level set methods have also been applied to character recognition 
(the determination of true character boundaries) and image enhancement in the presence of 
noise. One of the principle strengths of level-set methods is their ability to handle difficult 
geometries including sharp corners, cusps, and topological changes [13]. 

In the general level-set based shape generation case, a domain of possible shape geome- 
tries Q. 3 C R d in which the entire family of shape boundaries will lie is introduced where 
d (= 2 or 3) is the dimension of the shapes to be generated. For the generation of two- 
dimensional geometries ( d = 2), a generating function £7(x,tu) : —¥ IR parameterized by 
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Figure 2-4: Mesh plots of the level-set generating function Q{x.y\u) for u = 0.85 (left) 
and u = 0.60 (right). The resulting, closed shape geometries are shown shaded below each 
generating function. 


the vector = (i . . . , ) G -4 is selected where the Cartesian coordinate directions are 

x = ( z\ y ), the shape generation domain is Q 6 C 1R 2 , and the input parameter domain is 
A C The boundary of the shape geometry. is defined by the level-set 

= {x E | <?(x; u)) < 0}, (2.73) 

and is parameterized by the vector a In words, the body B ^ consists of all points in the 
x-plane for which the level-set generating function Q(x;uj) is less than zero. For a properly 
chosen generating function, this set is closed over the range of uj G A although it can be 
disconnected. 

To illustrate the level-set shape description described above, the sample generating 
function 

£(x; w) = u/i - e-K*- 1 ) 2 -*-* 2 ] - e -[ (x+1)J+ ^, (2.74) 

defined over the shape domain = [—2.0, 2.0] x [—1.0, 1.0], and parameterized by the 
single variable uj = is selected. The parameter uj\ serves to select particular geometries 
from the family defined by Equation 2.74 for the full parameter range uj 6 A = [0.5, 1.0]. 
Mesh plots of the generating function and shaded, two-dimensional plots of the resulting 
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shape geometries are given in Figure 2-4. The single shape parameter is set to uq = 0.85 
for the body on the left in Figure 2-4 and uq = 0.60 for the body on the right. For even this 
very simple case, the connectivity of the resulting geometry B T can change discontinuously 
for a continuously changing shape parameter u>. This is a critical result for the ultimate 
application of surrogate-based shape optimization explored later. 

The very simple level-set shape generation procedure described above can be easily 
extended by introducing a more general generating function. One approach is to treat the 
general form of the exponential function on which Equation 2.74 is based 

0i(x; u) = /3 ie -U*- x “) 3+ (*-*‘) 2 l (2.75) 

as the core for a more complex level-set function that is results from the application of 
superposition techniques. The level— set functions formed in this manner describe more 
general shapes. Each base exponential function has an associated center ( x ci ,y ci ) and 
strength 3 t . For N e exponential functions, this yields a generating function of the form 

S(x:u>) = fa - 5>(x:u,) = 8 0 (2.76) 

1 = 1 1=1 

where the vector of shape parameters is u> = {(/3 t .x rj . y Cj ) \ i = 0, . . . . N e .j = 1. .... iV t }. 
With the additional exponential functions and the greater flexibility afforded by the shape 
parameters, a much richer family of geometries can be defined. 

The level -set function in Equation 2.76 was explored for the eddy-promoter problem. 
It produces smooth bodies that can change topologies (from a single body to more than 
one) and also satisfies the preference that the description procedure include cylinders. The 
family of geometries that it produces, however, does not include more elongated, oval shaped 
bodies without increasing the number of exponential functions beyond what is feasible for 
accurate surrogate construction. For this reason, a different shape generation function, 
described in the next section, has been developed. 

2.5.2 Ideal Flow — Doublet Description 

Although not explicitly stated as such, level-set ideas have been used extensively in classical 
aerodynamic analyses [2]. The underlying assumption in these approaches is that the true 
flow solution over a body can be approximated by an ideal flow solution over the body 
[38] with an appropriate circulation superimposed to satisfy flow tangency at sharp trailing 
edges. The body surface in ideal flow analyses is a streamline of the ideal flow solution 
and general geometries can be analyzed by enforcing flow tangency on the body surface at 
discrete points[42]. 

Ideal fluid solutions can be expressed by a fluid potential function 4>i or stream function 
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t/ 7 , both of which satisfy Laplace's equation 


V 2 <f>i = 0, V 2 V>/ = 0. (2.77) 

The fluid velocities u = (u,u) are given by 


dcpi dipi 

(2.78) 

U dx dy 

drf>i dipi 

(2.79) 

dy dx 


The velocity solution u is for an incompressible, irrotational, inviscid fluid flow. Because of 
linearity, the superposition of flow solutions for which the velocities are given by 2.78 and 
2.79. and that that satisfy 2.77, will also satisfy the ideal flow governing equation. This 
superposition property is used to generate solutions over general body geometries. 


The ideal flow solution for the flow over a cylinder is described by the superposition 
of a uniform flow and a doublet. A uniform flow u(x.y) — U ^ is described by the stream 
function 


^uniform = ^ 


(2.80) 


A doublet is obtained by placing an equal strength source-sink pair a finite distance apart 
on the .r-axis and taking the limit as the distance between the source-sink pair goes to 
zero. The stream function for a doublet of strength A and center ( x c .y c ) is 


^doublet 


-M y - y c ) 

(x - x c ) 2 + (t / - y c ) 2 


(2.81) 


The stream function solution for the flow over a cylinder centered at the origin is obtained 
from the superposition of 2.80 and 2.81 for (x c .y c ) = (0.0). 


i cylinder /uniform , /doublet rr 

Wf =Wi +Wi = UocV — o— — 

x z + y z 


(2.82) 


The cylinder surface is defined by the t/y linder ~ 0 streamline (excluding the portions 
of the streamline on the x-axis) and the cylinder has radius /^cylinder = v/A/C/^. The 
relevance of this shape description to the original cylindrical eddy-promoter optimization 
problem studied previously [40. 77, 79] is readily apparent. By fixing the value of (7^, the 
doublet strength A then serves the sizing input analogous to d, the diameter of the cylinder. 
The position input, a, (measured from the lower wall of the channel) is the same as before 
giving the two-dimensional input vector Z^p nder = (A, a). 

The superposition idea described above can be further extended with the incorporation 
of a second doublet. The additional doublet allows for changes of topology (connectivity) 
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similar to the example level-set shape generation procedure described in 2.5.1 and the 
results of which are plotted in Figure 2-4. To accomplish this, one doublet of strength Ai 
is positioned at the origin and a second doublet of strength A 2 is positioned at (r d , 0). The 
uniform flow velocity is set to unity, U 0 c = 1. The stream function for this configuration is 
then 


M>l = y- 


Aiy 

x 2 + y 2 


My 


{x - r d ) 2 + y 2 ' 


(2.83) 


and the geometry is given by the ipi = 0 streamline with the segments on the z-axis 
removed. The parameters {r d . A^ A 2 ) determine the specific geometry that is selected. If 
the vector of inputs ( r d , Aj. A 2 ) € 9. DD is assumed to range over the input space 9 DD C JR 3 , 
then the connectivity of the resulting geometry (one or two distinct bodies) will depend on 
the values of (T- d ,Ai.A 2 ) and distinct regions of 9 DD can be identified that correspond to 
single-body geometries, fiP^, or to two-body geometries, fip®. and Q DD = 9.^ D U Q-P^. 

The original intent of defining a family of bodies that is a superset of the cylindrical 
configuration is satisfied for r d = 0. When r d — 0, a cylinder of diameter d. — \/A i + A 2 is 
obtained. The two-doublet generating function in Equation 2.83 is the form used for the 
eddy-promoter shape description described in the next section. 


2.5.3 Eddy-Promoter Geometry Generation Procedure 

The full set of geometry inputs for the two-doublet shape description are shown in Figure 
2-5. These include the strengths of the two doublets. A i and A 2 . the position of the second 
doublet center relative to the first, (x d .y d ). and the position of the pair in the channel 
measured from the lower wall. V . The distance between the two doublets in Equation 2.83 
is r d = Jx 2 + yl- The procedure implemented to obtain the geometrically meaningful set 
of inputs Z ep = {xd-.Vd-Y.X i,A 2 ) from the normalized, eddy-promoter, geometric design 
vector p z - {pz\<- ■ ■ ,PZb) G 9t z = [0. l] s is described below. 

The generation procedure defines the geometric design space f>z in terms of a set of 
geometrically meaningful limits put on the resulting eddy-promoter configurations. Because 
of the difficulty of solving problems on highly stretched finite element meshes, the eddy- 
promoter bodies are never allowed within <o of the upper or lower wall. The range of the 
of geometric input x d is set such that z m i n < x d < z max and a limit is set on the total 
doublet strength A r = Ai + A 2 of A min < A r < A max . Upper and lower limits are also set on 
the strengths of each individual doublet Ai min < Ai < A lmjn and A 2min < A 2 < A 2m ; n . The 
generation procedure that satisfies the geometric requirements listed above consists of the 
following 7 steps: 

1. The value of x d is set based on the normalized input p Zl \ x d - ^min ~^PZ l l-^max ^min)- 

2. The total doublet strength is set based on the normalized input pz 5 : A T = A min + 
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Figure 2-5: General geometry for the two-body eddy-promoter configuration. 


PZo(^max ^min) 

3. The total doublet strength distribution parameter, a, is set such that a m j n < a < 

Qmax* The limits are dependent on the total strength and are a m i n = and 

»max = l-«min* The distribution is set according to pz 3 : a = a min +pz 3 (a m ax-amin)- 

4. The strength of each doublet is set based on a and Xj\ Ai = aX? and A 2 = (1 — a)Aj. 

5. Set the value y d for the y-position of generator two relative to generator one based on 
Pz 2 : Vd = !/min+PZ 2 (ymax -J/min)- The limits on y d are dependent on the values of Ai, 
A 2 , and to and are y m \ n = -2(1 - i 0 ) + y/X\ + n/A *2 and y max = 2(1 - to) - >/Ai - n/X?. 

6. Generate the bodies based on the input values {x dr y d ). Ai , and A 2 and find the min- 
imum and maximum extent of the configuration in the y-direction. This requires 
placing a doublet of strength Ai at the origin, a second doublet of strength A 2 at 
(0 ? rd), finding the body surface, and rotating the configuration by 5. The maximum 
and minimum vertical extent of the eddy-promoter configuration, y^ in and yj^ ax . are 
used to set the final limits on Y , 

7. Set Y according to pz\\ Y — ^min +p^ 4 (y max + y m j n ). The limits on Y are functions 
of Vmm- 2/max: and l o and are r min = t 0 - y^ in and Y m ax = (2.0 - t 0 ) - y® ax . 

This somewhat involved set of steps must be performed to generate the bodies because 
simply positioning one doublet at (0,0) and a second at {x d ,y d ) will not produce a closed 
contour (or two closed contours) unless y d = 0. This requires the superposition of two 
doublets on the same axis and then the rotation to define the bodies in the channel. It is 
important to note that for the cases in which two distinct bodies are described, the resulting 
bodies are not circular cylinders. However, in the limit as the distance between the doublet 


46 



centers approach infinity, the bodies approach circular cylinders of radius R\ = \/\\ and 
/?2 — \/% 2 - 

Matlab [47] has served extensively for the shape generation. The steps described above 
have been coded into a Matlab routine that takes as input a sequence of normalized inputs 
p z . and writes as output a file that contains the surface node coordinates for each input, 
design point p The node coordinate file is read by the mesh generator to build the finite 
element solution mesh. 

As already mentioned, the primary benefit to defining the shapes in this manner is 
that it allows for topology changes relative to smoothly changing input parameters. This 
capability permits optimization over a wider range of configurations without having to resort 
to discrete optimization procedures. The most important characteristic of geometrically 
describing a shape space with more than one distinct topology in a level-set based approach 
such as this is that the input-output relationship can be expected to be continuous (and 
possibly smooth). This property is demonstrated empirically in Section 2.5.4. 

The level-set method described above can be extended through the incorporation of 
additional doublets to allow the definition of more geometrically rich body profiles. Ap- 
proaches such as this are very common in ideal flow analysis of airfoils. Finally, through the 
careful selection of generating functions that consistently define bodies for arbitrary (x t .y t ). 
the method could be extended to still more complex topologies. 

The above described steps generate a set of shapes that obey the geometric constraints. 
It also ensures that that the entire normalized input space 9-z corresponds to feasible 
designs. In the next section, an example that illustrates the topology change is given. The 
example shows the smooth transition from one to two distinct eddy-promoter bodies, and 
gives empirical evidence that the input-output functions are continuous across a topology 
change. 


2.5.4 Input-Output Response Example 

To demonstrate, at least empirically, that the relationship between the normalized geometric 
inputs, pz ■ and the outputs is continuous across a topology change, the finite element 
simulation described in Chapter 3 has been used to approximate two of the outputs. Only 
one of the normalized inputs, pz i, of the full normalized input vector p z has been has been 
varied. The first component of the normalized geometric input vector corresponds to the 
physical coordinate Xd shown in Figure 2-5. The bounds on x d have been set to x min = 0.00 
and x max = 0.75 and pz\ € [0. 1]. A total of 21 geometries have been generated with the 
values of pz\ clustered near the value at which the topology change occurs. The remaining, 
normalized geometric inputs have been fixed such that pz = (pz\- 0-60. 0.50, 0.50. 0.50). The 
lower and upper limits of the geometric quantities that correspond to fixed normalized inputs 
are A min = 0.109. A max = 0.109. a min = 0.85, and a max = 0.85. With the above information. 
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the geometry generation procedure described in Section 2.5.3 has been used to compute 
the surface coordinates for each input pz Representative eddy-promoter configurations 
produced by this set of inputs are plotted in the top half of Figure 2-6. 

For the computational case examined here, the Reynolds number has been set to Re = 
250. The two outputs of interest are obtained from the finite element simulation and are 


6 


FE 

0 


= tf E (Pz) 


1 

Nu(Zep{Pz), Re = 250) ’ 


(2.84) 


4>q E = iPo E {Pz) =< /(<; Zep(Pz), Re = 250) > (2.85) 


where Nu is the Nusselt number given in Equation 2.51, and f(t;Zpp, Re = 250) is the 
time-dependent forcing term in Equation 2.25 required to achieve the prescribed flow rate. 
The results of the computational study are given in Figure 2-6. Selected eddy-promoter 
geometries from the sample set are plotted in the upper half of the figure. The output 
responses versus geometric input value xj are plotted in the lower half of the figure. The 
normalized input pz\ is a simple linear scaling of x^, therefore, the output response plots 
versus pz\ would tell the identical story to those in Figure 2-6. In the output response 
plots, the dashed line indicates the approximate value of x d for which the topology change 
occurs. Single-body configurations have Xd values to the left of the dashed line and two- 
body configurations lie on the right. The results plotted in Figure 2-6 support the assertion 
that the input-output relationship remains continuous across topology changes. In fact, the 
first derivative appears to be continuous as well, further supporting (empirically) that the 
geometric description used for the eddy-promoter problem is well suited for the surrogate- 
based optimization techniques applied in Chapters 4 and 5. 


2.6 Reduced Problem for the Eddy-Promoter Example 

In this section, the particular set of inputs used for the eddy-promoter example problem 
and the corresponding design space are described. Second, the outputs of interest and 
the inputs of which they are functions are restated. Finally, the performance metrics that 
characterize the design preferences for the eddy-promoter design problem are presented. 

2.6.1 Inputs and Input Domain 

The full set of physical inputs are the variables that specify the specific geometry of the 
eddy-promoter inclusions Z ep — {xd-.Vd,Y, Ai,A 2 ) as well as the nondimensional, inverse 
channel height r/£. For the eddy-promoter problem examined in the Chapters 4 and 5, 
a reduced set of inputs is used. The generator strengths, A x and A 2 , are fixed and the 
remainder of the geometric variables, xj, yd, and Y serve as inputs. This gives a reduced 
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Figure 2-6: Input-output response results for a topology change. Selected eddy-promoter 
configurations are given in (a)— (f) and the output responses are plotted in (g) and (h) 
versus the geometric input x d . 
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set for the normalized input vector 


p = (pi,...,Pi) € Q = [0. 1] 4 (2.86) 

which has M = 4 inputs and is a member of the normalized design space f2. The fixed 
generator strengths are Ai = 0.0871, and A 2 = 0.1065. The relationships between the 
normalized inputs p and the geometric quantities (zrf. Y, r)i) are 

x d = ^min 4" pi (^max “ ^min)? (2.87 ) 

Vd ~ Vm’in 4“ P2(pmax ~ ymin): (2.88) 

^ — ^rnin 4“ P3(^max ~ ^min): (2.89) 

VL VLmin d“p4(pLmax — 'HL min)^ (2.90) 

where *£niin — 0.0, *^max — 3.o, 2/rnin? Pma x? ^rnin? and &re each determined based 

on to = 0.15 and the generation steps described in Section 2.5.3. 7lL m \n = 13.332, and 
77 /^max = 106.656. The limits on the inverse channel height correspond to a range, for 
the full heat exchanger, from 4 periodicity lengths. L through 32 periodicity lengths. The 
periodicity length is held fixed at l = 6.666. 

For convenience in the latter sections of the thesis, it will be necessary to make a 
distinction between the inputs for which the performance metric response must be modeled, 
either by the finite element simulation or by surrogates, and the inputs for which the 
performance metric response is known analytically. The vector of modeled inputs p m 6 
= [0. I] 3 consists of the first three components given in Equation 2.86 

Pm = iPl-.P2.Pi)-, (2.91) 

where the relationships in Equations 2.87 — 2.89 and the steps described in Section 2.5.3 
are applied to relate the normalized inputs to the geometric quantities. The analytic input 
p a E = [0, 1] corresponds to the inverse channel height 


Pa = (P4). (2.92) 

where is the identical input to that given in Equation 2.86 and which can be related to 
7fi through Equation 2.90. 

The design domain is likewise decomposed into modeled, Q m = [0, l] 3/m , and analytic. 
fl a = [0, l] 3/a , subspaces of fi where the number of modeled inputs is M m = 3 and the 
number of analytic inputs is M a — 1. The full design domain is the tensor product given 
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by 

fi = Q m x S2 a , (2.93) 

and the full design vector is 

p = (p m ,p a ) G ft. (2.94) 

2.6.2 Outputs 

The two-body, eddy-promoter problem has two outputs of interest. The outputs are each 
functions of the modeled design inputs, p m . The first output is the time- average of the 
nondimensional pressure gradient given in Equation 2.40 and expressed as 

tMPm) =< F(Pm.t) > = < fit) > (2.95) 

where as before. < • > refers to a time-average over a sufficiently long time. 

The second output is the inverse Nusselt 

0 O (P rn) = AT(P m) - v l ~ 7 (2 96) 

Nu(p rn ) 

and is defined in Equation (2.51). 

2.6.3 Performance Metrics 

The performance metrics have been described in detail in Section 2.3 and will only be 
restated here. The metrics have been selected to characterize the design preference for the 
full eddy-promoter heat exchanger, a segment of which is shown at the top in Figure 2-2. 
The engineering goals of the design problem are to transfer as much heat as possible into 
the channel and to do so with as little pressure drop as necessary. 

The performance metrics are functions of the modeled inputs p m through the output 
functions ^(p m ) and #o(Pm) and. although not the case for this problem, could be explicit 
functions of the modeled inputs as well. The performance metrics are functions of the 
analytic inputs p a explicitly. In fact, the analytic relationship between the performance 
metrics and p a is known. 

The first performance metric is the nondimensional pumping power 

'P(p) = ^(^(PmlT) = log 10 [i/>o(P (2.97) 

which is a function of the output </> 0 (p m ) and of the analytic input p Q through the nondi- 
mensional. inverse channel length r /£. The response of 'P(p) to the input p a is known 
analytically through the normalization relationship in Equation 2.90. 
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The second performance metric is the nondimensional wall temperature difference 


©(P) = ©(0o(Pm),P) = lo glO 


1 n , . 1 

— — — - + 0o(Pm) — 
RePr T) L J 


(2.98) 


which, similarly to ^(p), is a function of the output #o(Pm) and of the analytic input p Q 
through the nondimensional, inverse channel length r)i. The response of 0(p) to p a is also 
known analytically through Equation 2.90. 
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Chapter 3 


Numerical Approach 


In this section, the numerical approach used to solve the governing equations 2.25, 2.26. 
2.31 is described. A second-order accurate in time, semi implicit, split, time advancement 
scheme has been used [23]. The spatial discretization uses a Galerkin finite element approach 
with triangular elements. Second order, isoparametric elements are used to approximate the 
velocity and temperature solutions and first order, linear elements are used to approximate 
the pressure. A variable time-stepping scheme has been employed to allow the time-step 
size to vary as the solution is advanced so that the Courant condition is satisfied. 

In the following sections, the details of the numerical approach are presented. First, the 
spatial discretization is described including the elemental matrix treatment, the strategy 
used to form the convection operators on-the-fly. and the required quadrature accuracy to 
ensure that heat does not build up. Second, the second-order accurate, temporal discretiza- 
tion is presented, the method used to achieve the prescribed, constant flowrate is given, and 
the iterative and direct solvers used to invert the matrices are described. Finally, a sample 
of the code performance for a representative problem and for the various solver options is 
presented. 

3.1 Finite-Element, Spatial Discretization 

For symbolic expediency, in this section the velocity vector components will be expressed 
indicially, u = (ui,U2). This runs counter to the (u, v) form used in the remainder of 
this thesis but greatly simplifies the presentation of the numerical implementation. The 
variational weak form of the governing equations 2.25. 2.26, 2.31 over the solution domain 
Q for the solution variables u(x, t) £ [//i 0 (fl)] 2 . p(x. t) £ T" t0 (n). and 9(x.t.) £ //.L(fi) is 

(t/i. + (vi.u • Vtq) + Vu,-) - (V^p) = (vi.fi), i = 1,2. 

Vv £ [//i 0 (f>)]*. (3.1) 
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(V • u. g) = 0, v ? el 2 # o(«) 


(3-2) 


and 

(w,0<) + (i«:U’ (V0)) + ^ e p r ( Vu; ’ v<? ) = -7(“Mi-et) + fep r ^ dn ‘( w )' 

Vw € H l # 0 (fi), (3.3) 

where 

= [ wdl, (3.4) 

./an, 

Q, is the lower wall of the computational domain. L^ 0 (f?) is the set of all sq-periodic 
functions that are square integrable, and is the set of all aq-periodic functions that 

are square integrable and have first derivatives that are also square integrable. The product 
(a, b) is defined as 

(a. b) = [ abdx . (3.5) 

Jo. 

A finite element, spatial discretization of the governing equations is pursued. A IPi-IPx 
(Taylor-Hood [70]) elemental basis is selected. With this choice, the velocity and tempera- 
ture solutions are represented by second order triangular element and the pressure is treated 
with linear triangular elements. The reason for choosing the Taylor-Hood elements is that 
it is known to satisfy the inf-sup condition [33] for the existence of a unique solution and 
can be conveniently implemented. The second order node points consist of the vertices of 
each element as well as the midpoints along each element side and the linear ’’pressure'* 
mesh is simply a subset of the second order mesh consisting of only the node points at the 
vertex of each element. 

The full solution domain Q is subdivided into M el triangular elements such that 

X el 

n= (J (3.6) 

k e = 1 

and is the triangular element with index k € . The discrete solution to the governing 
equations 3. 1-3.3 can then be expressed as: Find u h (x.t) G [V^] 2 , p h (x,t) € P^, and 
Q h (x,t) G Vh where the following are satisfied 

f)~.h i 

(v i ^) + (v i ,u h -Vu*) + —(Vv i ,Vu*)-(Vv l .p h ) = (v l jh, * = 1)2, Vv6[^] 2 , (3.7) 

(V-u\g) = 0, V <7 € (3.8) 

+ (w,u h • (V0*)) + jX-(V w .V9 h ) = —y(iv. u h -ei) + wdS , 

Vw e V h . (3.9) 
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The superscript l h' refers to the discrete solution and the discrete subspaces are defined as 

p h = Wn*. 6Fi(fi*‘)}ncJ 0 (fi), (3.10) 

v h = {t>| n * e e F 2 (fi fc ‘)}n/fJ 0 (n), (3.u) 

where C# q (Q) is the space of all xi-periodic functions that are continuous and have zero 
average. 

A Galerkin finite element approach, in which the solution and test function bases are 
identical, is pursued. A global basis for the velocity and temperature is introduced 

&(xj) = Sij, 1 < i,j < V h = span{<&. i = L . . . ,Af v }, (3.12) 

where jV v is the number of global second order nodes and Xj is the global coordinate of 

node j. The basis is used to express a general function v over the full solution domain fi as 

A' 1 ' 

v b {x) = v n 4> n {x) (3.13) 

77= i 

where v t = v(x t ). Similarly, a second basis for the pressure is introduced 

ipi (Xj ) = Sij, 1 < i.j < M p . Ph = span {01- i = 1 ? A/ rp }. (3.14) 

where Af p is the number of global pressure nodes. The basis is used to express the general 

function r/. discretely over Q as 

A' p 

q h {x) = ^< 7 «tMx) (3.15) 

n = 1 

and. as before. q t = q{x t ). 

To proceed with the discretization, the solution variables and the test functions are 
expanded with the appropriate basis function defined above. The expanded forms 

,V v A' 1 ' 

U h = ]T U n 0 n (x). V = V„0 n (x). 

n=l n=l 

A* p A* p 

P k = Y^Pn^nix), 9=^9n^n(x). (3.16) 

n=l n= 1 

A* 1 ' A' 1 ' 

= ]T0 n <Mx), w = Y^ u , «^(x). 

77 = 1 n=l 

are substituted into Equations 3. 7-3. 9. The resulting spatially-discrete. governing equations 
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are obtained 


(3.17) 

(3.18) 


B-^l + C(u)u, + — At i, - D Jp = B t = 1.2. 


D,u, = 0. 


B0 t + C(u)0 + — ^-A9 — — 7 B(u • e t ). (3.19) 

Kerr 

where, for convenience, the superscript 'h' has been dropped. The finite element matrices 
are 



(3.20) 

(3.21) 


D 


Uj - 




uV4>i4>j dx. 


(3.22) 

(3.23) 


The finite elemental matrices given in Equations 3.20 — 3.23 are evaluated elementally, and 
a direct stiffness summation procedure is used to apply the corresponding global matrix 
operators when needed. 


3.1.1 Elemental Matrices - 1P> Isoparametric Formulation 

At the elemental level, each element is transformed from the Cartesian coordinate system to 
an elemental coordinate system through a a transformation. For subparametric elements, 
an affine transformation is used. For such a transformation, the elemental coordinates 
£ = (s ; ■ s 2 ■ S 3 ) take the form 


£1 = ai + bis + c\y, 


£2 -0,2+ b>x + c 2 y. 


£3 = «3 + b 3 x + Ciy. 


where the coefficients are given by 


(x 2 y3 - *31/2 ) 

6l = 

(y 2 - y. 3 ) 


O 

II 

IO 

2 A ’ 

Cl = 

(^yi - X 12 / 3 ) 

b 2 = 

(y3 - yi) 


= 2.4 

2.4 ! 

c 2 = 

{xij/2 - x 2 yi) 

b 3 = 

(yi - IJ 2 ) 


“ 3 = 2.4 ’ 

2.4 ' 

C 3 = 


(x;j — x 2 ) 
2 A 

foi ~ ^3) 
2.4 

(x 2 — X]) 
2 A 


(3.24) 


(3.25) 


and .4 is the area of the element. Note that only two of the affine coordinates are independent 
(as expected for M 2 ) and the the relationship that si +£2 + £3 — 1 holds. If subparametric 
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Figure 3-1: Affine mapping from the x-plane to barycentric coordinates £ = (6-f2-f3)- 


elements are used, the expressions in Equations 3.24 and 3.25 allow the elemental matrices 
to be evaluated analytically. Because isoparametric elements have been used in this work, 
the elemental matrices can not be evaluated analytically and the relationship in Equations 
3.24 and 3.25 can not be applied. Instead Gaussian quadrature must be used to evaluate 
the elemental matrices' integrals. 

Over element k e . a function can be expressed in terms of the elemental, linear basis 
functions as 

«*■«) = Erf- MCI. o- 26 ) 

1=1 

where the linear basis functions are 

/ll(f) = fli ^>(0 = 6i ^3(f) = s3- (3.27) 

and qf e is the function value at node i of the element k e . Similarly, a function can be 
expressed in terms of the elemental, second order basis functions as 

= (3-28) 

;=i 

where the second order basis functions are 

Mf) = f i(2fi-i), Mf) = 466- 

MO = 6(26 - 1)- Mf) = 466- (3-29) 

/ 13 (f) = 6(26 - i), M 0 = 466- 

The derivatives of the basis functions with respect to the barycentric coordinates. 6 are 
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obtained from the expressions in Equation 3.29, and the relationship that £3 = 1 — £1 + 6 
and are 


dhiiO 

56 

gMO 

56 

gMg) 

56 

5M6 

56 

5M6 

50 

gMO 

50 


46 - 1, 

5M€) n 
56 

46 - 1, 

5M0 n 
56 

4(6 +6) -3, 


46, 


-46, 

^ = 4(1-6-26), 

4(1-26-6), 

5M6 

0(2 = 4 '-" 


(3.30) 


A direct stiffness summation procedure is used to form the matrices given in Equations 
3.20-3.23 based on elemental components. This is accomplished by using the local basis 
functions given in Equations 3.27 and 3.29 and transforming the integrals to £-space. For 
example, the elemental contribution for stiffness matrix is 


/ 5M€)5Mfl , 5MO dhj(Z)_\ 

V dx dx dy dy ) 

where the Jacobian J(£) is given as 



(3.31) 


J(t) = 


dx dy 


d£\ d£ 2 


dx dy 

d £ 2 d£i 


(3.32) 


and is a function of the coordinate £ because of the isoparametric implementation. The 
(x,y )~ coordinates are expanded in terms of the basis functions in Equation 3.29 as 


6 

2=1 


6 


y(0 = 'jT,yihi(£). 
2 = 1 


(3.33) 


Plugging the coordinate expressions into Equation 3.32 gives 



which is now expressed entirely in terms of the £ coordinates. 



(3.34) 


The derivatives within the parentheses in 3.31 (of the elemental basis functions with 
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respect to x and y) are expanded as 


dHi) _ dhi{t)dix + dhiW db 


dx 


d£i dx d £ 2 dx ’ 


dhi(t) dhdS) dh + dhi(Z) d& 


(3.35) 


(3.36) 


dy d£ i dy d& dy 

To obtain derivatives of the elemental coordinates with respect to the Cartesian coordinates, 
the relationships 

dx = +i£ 2 d£ 2: (3.37) 

dy = 1 + y^efo, (3.38) 

are solved for the differentials d£\ and 2 . This gives the expressions 


dtx = ^ C '^ X — (3.39) 

_ ^£2 2/£l 

= xg. rf j- r (3.40) 

from which the desired results are easily obtained. Recognizing that the denominators in 
Equations 3.39 and 3.40 are the Jacobian (3.34). the expressions for the derivatives of the 
elemental coordinates with respect to the Cartesian coordinates are 


d£i _ fe 
dx j(ty 
d& _ _Vix_ 
dx J{ZY 


dji _ x ^ 
dy J[tV 
dj 2 _ 
dy J(€)' 


(3.41) 


Again, the derivatives of the Cartesian coordinates with respect to the elemental coordinates 
can be obtained, as they were for the Jacobian in Equation 3.34. by using the expressions 
in Equation 3.33. 


An expression for the stiffness matrix that is entirely in terms of the elemental coordi- 
nates. f . can be obtained by first substituting Equation 3.41 into Equations 3.35 and 3.36 

t00btam - ' gMfl dy dhYZ) dy \ 

a? 1 di 2 dt 2 dtJ' 

f dhi{£) dx 


dhi(Z) 


dx 

dhi{i) 


dy 


1 

m 

1 

m 


(3.42) 


dhi(£) dx 


)• 


(3.43) 


di 2 d^ d£ 1 d{ 2 . 

Finally, substituting Equations 3.42 and 3.43 into Equation 3.31 gives an expression for the 
stiffness matrix that is a function of £ only 


A- 


( 0 = f 

Jn k ‘ 


1 


n*. J(t) 


f dhj(£) dy 
V d£ i d£ 2 


dhj(£) dy\ ( dhyt) dy 

dti 2 dh ) V dt-2 


dhjij) dy \ 

di 2 dtJ 
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dt (3.44) 


/3/i;(4) dx dhi(£) dx \ /dhj(£) dx dhj(£) dx V 

V d& d£i d£i d£ 2 ) V d( 2 d£i d&J _ 


where the Jacobian, J(£)> is given in Equation 3.34. 

The remainder of the elemental matrices are obtained by following the same procedure 
as was used for the stiffness matrix. The resulting expression for the mass matrix is 

B?/(0= / JMhiMhite) <%. (3.45) 


The gradient matrix in Equation 3.22 is split into two components; the first for the derivative 
in the x-direction and the second for the derivative in the ^/-direction. The resulting two 
matrices are 


L 


(dhi(i) dy c>MO dy 


\ 6»s C l 

dh k (£) dx 


d& d£i 

dh k {£) dx 


d& d‘i d& 


hj(Z) dt 

hit) 


(3.46) 


(3.47) 


and it is interesting to note that the Jacobians completely cancel in the above expressions. 


The elemental component of the full convection operator in Equation 3.23 is broken 
into its directional components. The velocity vector that is embedded in the integral is 
expanded with the local basis functions similarly to Equation 3.33. The expansion allows 
the velocity values at the nodes to be taken out of the integral. The elemental convection 
operator is 

C-; <u*‘) - £ (ul)-' c4 + £ (u 2 tc 2 *j k , (3.48) 

1=1 1=1 

where the convection operators take the form 


C 2 &(€) 


[ hiMhjd) 

Jn k ‘ 

[ hi{t)hj(t) ( 


dh k (£) dy 
<9£i d£ 2 

dh k (£) dx 

d& d£ i 


dhkiX) dy 
di 2 d£i 



dhkjj) dx 
d£i d£ 2 



(3.49) 

(3.50) 


For subparametric elements in which the second order nodes are at the midpoints of 
the vertices, the elemental matrices in Equations 3.44-3.50 can be evaluated analytically. 
However, subparametric elements introduce ; skim [69] errors along curved boundaries. To 
reduce the skin errors, isoparametric elements are used along the curved boundaries. With 
the isoparametric elements, it is not possible to compute the integrals of the elemental 
matrices analytically, and Gaussian quadrature is used instead. Using Gaussian quadrature, 
the integral of a function g(£) : Q ke — » IR is approximated as 



\ rq 


i=i 


(3.51) 
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where area** is the area of element k e . The weight at point i is s;* and the position 
to evaluate the integrand is quadrature point W hen used for the matrix evaluations, 
the area of the elements in £ -coordinates is always area* 1 * = 1/2. A 5 th ~order Gaussian 
quadrature scheme in which Af q = 7 has been used for all of the simulations. The weights. 
zu t , and quadrature points, are given in Appendix C.l. 


3.1.2 Convection Operator Treatment 

The elemental stiffness, mass, and gradient matrices in Equations 3.44-3.47 are computed 
for each element k e = 1. . . . , A f ei and stored in three-dimensional arrays during the initial- 
ization phase of each run of the Navier-Stokes code. Neglecting symmetries, the storage 
required for each of the stiffness matrix, mass matrix, and the two gradient matrices com- 
bined. is 6 x 6 x j\f el . which is not unreasonable from a memory standpoint.. However, each 
convection operator requires 6 x 6 x 6 x ;V ^ for storage which has a dominant impact on 
the code memory requirement. While this issue is only moderately significant for the two- 
dimensional solver, it dominates the storage for three-dimensional problems and severely 
limits the size of problemsthat can be solved. 

Fortunately, the application of the convection operator can be evaluated for each time- 
step as needed (on-the-fly) at the expense of only a moderate increase in computational 
effort. This completely eliminates the overwhelming storage penalty. To demonstrate the 
scheme used, the convection operator term in the u - momentum equation is used as an 
example. The velocity vector in two-dimensions is u = (u\« ^ 2 )* The variational weak form 
for this term is 


(u\ u- V-itO 



(3.52) 


where w is the test function. The test function and the velocities are expanded in terms of 
the elemental basis functions 


w 




(€) = XX 

;= 1 


t*? e (o = y>i)f'M$). «2*(o = i> 2 )f e MO- 


(3.53) 


1=1 


1=1 


where iu^ e , (wi)f'- and (iti)** are the values of the test function and the velocities at node 
i of element k e . 

If only the first term on the last line of Equation 3.52 is considered first, the expansions 
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in Equation 3.53 and the result given in Equation 3.49 can be used to obtain 


6 6 6 

L = c 4«>- 

j9 - k ' 0x i=\j=\k=\ 


(3.54) 


Similarly, using the same expansions and the convection operator in Equation 3.50, the last 
term on the last line of Equation 3.52 can be written as 


6 6 6 


(3.55) 


L j{ z )WU2 i£ d * = EEWH-(“i)!' c 4^ 

ju * uy t=ij=ijt=i 

Recognizing the the convection operators in Equations 3.49 and 3.50 are of the form 

C ■&«) = C4({) = f^' GjjKS) d(. (3.56) 

'dh k (Z) dy dh k { £) dy 


where 


G%{Z)=h l {i)h j {S) 


G%{Z)=h l (Z)h j (Z) 


i d& d^2 d£i 
dhkiO dh k (£) dx 


(O fa \ 


(3.57) 


(3.58) 


0*2 d£i d£i 

and evaluating the integrals by Gaussian quadrature as in Equation 3.51 gives the following 
expressions 


1 


,\' q 6 6 6 




/-I 1=1 J=\ k=l 
6 6 6 


L d s* l E £ E £ ^“'.■ e («i)j , («2)^GS 2 *(€/). 

• /u e y ^ /=l 1=1 j=U-=i 


(3.59) 


(3.60) 


where is the coordinate of quadrature point l. 


The summations over i in Equations 3.59 and 3.60 are dropped because the solution 
will be for any test function inf' € V/,. The expressions for the convection operators could 
be implemented in a code with nested summations exactly as suggested by Equations 3.59 
and 3.60, but this would be extremely inefficient. Instead, by expanding the expressions for 

/-I r* 

G t ;j' k {£i) and G ( d.(^). the terms in Equations 3.59 and 3.60 precomputed and stored inside 
of the direct stiffness summation procedure. This greatly reduces the number of operations 
from a brute-force implementation. 


First, for Equation 3.59, the regrouping takes the form 

V"(j g g 

/= 1 j= 1 k=l 
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.(3.61) 



At this point, all that remains is to clean up the terms in Equations 3.61 and 3.62. 
The j- and ^-summation terms in Equations 3.61 and 3.62 can be precomputed at each 
quadrature point and stored. The terms are 


6 





(3.63) 



6 

J=l 

(3.64) 


= Z( n i)*' 

f0hj((,)0y dhj(£i) dy \ 

(3.65) 


ii 

Isj 

>r 

(dhj{£,)dx dhj(£i) dx \ 

V dt 2 dti dti db) ‘ 

(3.66) 


where = ((uOf*,... ,(«i)J 4 ) and = ((u 2 )f e -- • • . («2)£ e ) are vectors of the velocity 
components at each node of element k e . A vector .N' n long is required to store each of the 
terms in Equations 3.63-3.66. Additionally, the values of the basis functions and the deriva- 
tives of the basis functions can be precomputed at each quadrature point, further reducing 
the operations needed during the time advancement. The only elementally-dependent quan- 
tities that are needed for the above calculations are for U*', u k ', and the derivatives of x 
and y with respect to the elemental coordinates. The derivatives are computed as they are 
for the expression in Equation 3.34. 

Finally, inside of the loop over the elements of the direct stiffness summation procedure, 
the terms / c (u* e ,€/), I c {v l£{v k ',Z t ).. and Jf(ii*‘,£t) are 
computed and stored in A/’ 9 -length vectors. For the elemental computations, the terms 
I c (u ke ,£i) and I c {v k ',£i) each require 2 x 6 x A f q operations and the terms (u ke ,€i). 
I^(v k ‘.^ t ). I C (u ke ,£[), and Iy(v k€ -£i) each require 5x6 x jV q operations. The derivatives 
of the spatial coordinates with respect to the ^-coordinates in Equations 3.65 and 3.66 
require 8 x 6 x M q operations. The effect of the convection operator is added to the the 
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right hand side of the discrete equations with a direct stiffness summation procedure. The 
elemental contributions to the right hand sides of the x- and y-momentum equations are 


1 A ' q 

RHS !*' = + , (3.67) 


l— l 


1 

RHS^ = + , (3.68) 

z Z=1 

respectively. 

The total operation count required to apply the stored convection operators to both 
momentum equations is 1512 x,\f el . The total number of operations required to compute the 
right hand side contributions in Equations 3.67 and 3.68, and to apply them globally to both 
momentum equations using the direct stiffness summation procedure, is 2712xN q Af el . For a 
5 th -order quadrature scheme in which A f q = 7, the on-the-fly evaluation requires roughly 13 
times as many operations to apply convection contribution full-storage technique. However, 
the memory requirement is reduced from 2 x 6 3 xA re/ for the full-storage technique to roughly 
30 x j\f q for the on-the-fly evaluation, which is trivial compared to storage needed for the 
other elemental matrices. 


3.1.3 Temperature Solution Periodicity 

Care must be taken in the solution of the governing equation for the periodic part of the 
temperature given in Equation 2.31. For the continuous equation, it has been shown in 
Section 2.2 that the last term on the right hand side of Equation 2.31 will ensure that there 
will not be a build-up of heat over long periods of time. This is not guaranteed for the 
discrete case. If. however, a sufficiently accurate Gaussian quadrature scheme is used, it 
can be shown that the discrete equations will be exactly satisfied by the value of 7 given in 
Equation 2.39. 

For convenience, the weak form of the periodic temperature Equation 3.3 is repeated 

{w,e t ) + (w, u • (V0)) + -^T- (Vu;, Vfl) = —j(w. u • ei) + — ■ ■ — ^,(w). (3.69) 

A procedure similar to that used in Section 2.2 to find Y is followed for the discrete equations 
to determine the required accuracy for the Gaussian quadrature rule so tha^fbr 7 = ±p e p T , 
heat does build up over long time periods. Integrating the equation over the discrete solution 
domain Q is equivalent to setting the test function w = 1 . 

First, the requirement that heat not build up over long time periods is equivalent to 
requiring that 

1 C +r 

<(u;, 0 t )>=- / (u>, d t ) dt = 0 , (3.70) 

t j t 
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for w = 1, and r on the the order of the period of the steady-periodic state. In the discrete 
case, the time-average in Equation 3.70 is equivalent to summing the 2 nd -order, implicit, 
backward difference stencil over a sufficient number of time-steps. Because the coefficients 
of the stencil given in Equation 3.83 sum to zero, inner terms of the summation cancel 
exactly, and Equation 3.70 is satisfied discretely when averaged over enough time-steps. 


For the second term on the left hand side of 3.69 to be zero when w = 1, the quadrature 
rule must be sufficiently accurate so that the integral over each element for w = 1 is exactly 
evaluated. The test function w will be is included in the elemental integrals. At the end of 
the analysis, setting w = 1 will be equivalent to evaluating the integral over the element. 
The integral in terms of the elemental coordinates £ is expressed as 


[ icV • (u0) dx. = [ wJ(£)V ■ (u0) d£. 
Jn k ‘ 


(3.71) 


The solution values and test function inside of the integral are expanded according to the 
elemental basis functions in Equation 3.28 to obtain 




(3.72) 


(u0)(O = X>0)f'M£). 
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(3.73) 


(3.74) 


(3.75) 


~~dT"tx n Q y ' ,l m\ % % dbj 7 

Note that a modified expansion has been used in which the velocity-temperature prod- 
uct is treated as one variable. This form has been implemented in the code and allows a 
lower quadrature rule be used than if the velocity and temperature were expanded sepa- 
rately. By substituting Equations 3.72 — 3.75 into Equation 3.71, the convection term of the 
temperature equation is written as 


jv-(u9)dx = f E(^)i' e £>xI(4) + E( v0 )f e£> i'(O 

• /n * e j = 1 li=l i = 1 

6 6 

= £ w k / J>0)?' f hj(Z)D k x Kt) d£ + 

jtl Jnk ' 

6 6 

J hjMD^iOdt, 


(3.76) 


]= l i=l 
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where the substitutions 


£*'(0 = 


dhdQ dy c )hj(t) dy 

di\ d& d£x ’ 


D k e( c\ ... gMjl — dx 

yi di 2 dfr 


(3.77) 


(3.78) 


d£i 9^2 ’ 

have been made. For w = 1, the basis functions /ij(£) in the integrals of Equation 3.76 can 
be neglected. The highest remaining power, in terms of the components of £ = (£i,£ 2), in 
the integrals are of 2 nd -order. This implies that the integral will be evaluated exactly for 
each element, and subsequently for the summation over all elements of the domain, when 
a Gaussian quadrature scheme of at least order two is used. As stated in Appendix C.l, a 
5 th -order quadrature scheme has been used that satisfies this requirement. 

The third term on the left hand side of Equation 3.69 goes to zero for w = 1 from 
taking the gradient of the constant function, Vw = 0. All that remains is to show that the 
two terms on the right hand side exactly balance. The first term on the right hand side is 
integrated over the solution domain Q by setting w = 1. Evaluating the discrete integral 


~7(tc — L u ei) = —7 QL 


(3.79) 


where the result is dependent on the flowrate. The flowrate is prescribed Q = | and is 
exactly enforced as described in Section 3.2.1. 


The final term on the right hand side of 3.69 is evaluated for w = 1 to get 


1 


RePr 


?an,{w = 1 ) = 


l 


RePr 


(3.80) 


The sum of Equations 3.79 and 3.80 must be zero. Solving for 7 gives the desired result 

3 


7 = 


4 RePr 


(3.81) 


The accuracy for the Gaussian quadrature rule necessary to evaluate the discrete inte- 
grals exactly, and to ensure that for 7 given in Equation 3.81 that heat will not build up in 
the periodic domain is determined by the integrands in Equation 3.76. The quadrature rule 
must be such that a second order polynomial function is integrated exactly. A 5 th -order 
scheme has been used which meets (exceeds) this requirement. 


3.2 Temporal Discretization — Split Time-Stepping 

The spatially discrete, analytic in time, governing equations are given in Equations 3.17- 
3.19. The temporal discretization pursued is a semi-implicit, fractional time-stepping 
scheme [23, 28]. This approach reduces the time-advancement problem to a sequence of 
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smaller subproblems and allows the convection operator to be treated explicitly while the 
remainder (the Stokes subproblem) is advanced implicitly. In the scheme that has been 
used here, a 3 rd -order accurate, explicit Adams-Bashforth integration scheme is used to 
advance the nonlinear convection terms, and a 2 nd -order accurate, implicit backward dif- 
ference scheme is applied to the Stokes subproblem. 


A variable time-step formulation is used. The reason for this approach is that, due to 
the unsteady nature of the solutions, the stability criterion is likely to vary significantly 
as the solution is advanced. It is difficult to determine a priori a time-step that will be 
stable for the entire solution evolution, and it is inefficient (from a computational effort 
standpoint) to use the minimum step-size for the entire solution advancement even if it 
could be determined. The semi-implicit form for the velocity equations is 


B«,u" +1 + ( 2 u? + (3U"- 1 ) + Ta„" +1 - DlV" 1 


= £,?,C(u"-«)u"-» + Bf; 

< 7=0 


n-t-1 


i = 1.2. 


(3.82 


where the 2 n<l -order. backward Euler coefficients are 


2Af„ + A/„-i At. n + At. n -i _ A/,, 

At n ( At n -r- At„_ ] ) * A/„ At,, i At n {At n + At„-i) 


(3.83) 


As expected, for the constant time-step case where Af„_i = At n - At. the coefficients in 
Equation 3.83 reduce to give the well known finite-difference stencil for the first derivative 


g’(t) 


n + 1 


3 g n+1 -V +g n - 1 
2A t 


+ 0{At 2 ). 


(3.84) 


The coefficients fi q . q = 0. 1.2, have been derived in Appendix C.2 and are 
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(3.87) 


The terms in Equation 3.82 can be collected to obtain 


2 


Hu- +1 - D Jp n+1 = -B(C 2 u" 4- C3<- l ) - E /3 g C(u n -"K 

< 7=0 


n-q 


+ Bf n+1 . 1 = 1.2, (3.88) 
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where H is the second-order, Helmholtz operator 


H=^A + Cl B. 


(3.89) 


The fractional splitting scheme approximates semi-implicit form in Equation 3.82 with 
a sequence of smaller steps. The result is 

2 

Hu n+1 =D r p n+1 -B(C 2 u" + C 3 U"" 1 ) - ^/3 9 C(u"- ,? )u"- <? + Bf n+ \ (3.90) 

<?= 0 

Epcorr = -ClDu n+1 (3.91) 

U n+1 = u n+1 - -^-B _1 p corr . (3.92) 

Ci 

p n * 1 = jp l -hp corr (3.93) 

where u is an intermediate approximation to the velocity at step n-f 1. E is the consistent 
Poisson operator 

E = + D 2 B-W. " (3.94) 

Pcorr = p n ^ X — is the pressure correction term, and jP^ 1 is an approximation to the 

pressure at time step n -f- 1 and is simply a linear extrapolation from the previous two time 
steps 

p n - l =p n + (p n ^p n - 1 )-^. (3.95) 

The solution steps of Equations 3.90-3.93 are as follows; (3.90) approximate the pressure at 
n 4- 1 and solve for an intermediate velocity u ri ~ l , (3.91) solve for a pressure correction term 
Pcorr by enforcing divergence-free flow. (3.92) use the pressure correction term to correct 
the velocity approximation. (3.93) use the pressure correction term to update the pressure. 

The governing equation for the periodic part of the temperature, 9 , is decoupled from 
the flow equations, and is updated after the flow solution is updated. An identical scheme is 
used to advance Equation 3.19 in time. The resulting, semi-implicit, fully discrete equation 
is 

2 

H e e n+1 = -B(Q 2 d n + C 30 n_1 ) - Y. A 7 C(u n -‘ ? )0 n -9 - 7 B(u n+1 ■ e0, (3.96) 

7=0 

where H# is the temperature Helmholtz matrix 


H fl = 


1 

RePr 


A + CiB. 


(3.97) 


The inversion of the Helmholtz operator H in Equation 3.90 and of H# in Equation 3.96 
are elliptic solves and can be accomplished with an iterative, conjugate gradient algorithm, 
or with a direct, sparse LU-factorization. Because H and are well conditioned, Equa- 
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tions 3.90 and 3.96 can be solved very efficiently with a preconditioned conjugate gradient 
method. Inverting the E operator in Equation 3.90 is more difficult as the conditioning 
is very sensitive to the homogeneity of the finite element mesh. For solution schemes in 
which all inversions are iterative, it requires the bulk of the solution time. Because E is 
much smaller than the other matrices (its size is determined by M p < Af v ) and its inversion 
dominates the solution time, it is the best candidate for treatment with a direct solver. 
Details of both the iterative and direct solvers used in the code are given in Section 3.2.2. 

In both 3.91 (inside of the Poission operator) and 3.92, the mass matrix B must be in- 
verted. For a spectral element spatial discretization, this is a trivial matter as the elemental 
mass matrices are diagonal. However, the mass matrix for the P 2 , triangular elements is 
not diagonal, and therefore must be inverted with either an iterative solver or direct, sparse 
solver. To avoid the additional matrix inversions resulting from the use of the full mass 
matrix, the diagonally lumped mass matrix has been substituted. The elemental, diagonally 
lumped mass matrix B^' for element k e is 


B* e = 


B 


Trace(B^' 




i = 1. 


. 6 . 


(3.98) 


where | f> Av j is the area of element f 2*' computed from the full elemental mass matrix and is 


in*' I - 


i=i j=i 


(3.99) 


3.2.1 Flow-Rate Specification — Green’s Function 

In this section, the method used to determine the time-dependent, pressure forcing term 
p+i _ (y n +i.o) in Equation 3.90 is described. The value of the / n_rl must be set to 
achieve the desired average flow velocity V — |. This is accomplished by first solving in a 
preprocessing stage for a Green s function u*. The Green s function is obtained by solving 
the 'Stokes-like’ part of the spatially discrete Equations 3.17 and 3.18 for a unit pressure 
forcing term. Then, the full set of governing equations are advanced in time, neglecting 
the pressure forcing in Equation 3.90, and, after each time step, superposition is used to 
achieve the correct flow rate and to determine the required value of / n+l . 

The preprocessing step (which only needs to be solved once) consists of solving the 
Stokes equation 

Hu* = Bl + Dfp*, i = 1,2. (3.100) 

for the velocity u* = (u^u^) and the pressure p*. The unit forcing 1 = (1.0) is a vector 
of ones on the z-equation nodes only. The finite element matrices are identical to those in 
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Equation 3.90. The flowrate is computed 


Q ' = hL u ' dx ' 


(3.101) 


and is used to determine the pressure forcing term f(t). The intermediate velocity u* and 
pressure p* solutions will be superimposed with the time-step solution to determine the 
correct Navier-Stokes solution. 

The second, time-advancement phase of the solution proceeds exactly as outlined in 
Equations 3.90-3.93, but with the forcing term neglected. The modified, intermediate 
solution phase has the steps. 


- B(C 2 u/" + Csu/"- 1 ) - £ 0 q C(u/"-*)u/ n -* 

(3.102) 

< 7=0 


E/W, = -ClDu7"+ l 

(3.103) 

U/ " +1 = ST +1 - 1b - l Pconr 

Cl 

(3.104) 

Pi n ~ l = Pi n + i +Pcorr, 

(3.105) 


after which, the pressure forcing is determined and the solutions are superimposed. The 
forcing term is simply 


r * 1 = 


Q 


n + 1 


Qi 


Q' 


where the intermediate flow rate Q/ is 


Ql = f u { dx, 

21 J-pEP 

and the flow rate at step n + 1 is the prescribed flow rate 


(3.106) 


(3.107) 


Q n+l = IV = If 
3 


(3.108) 


The final step is the superposition of the solutions to obtain the velocity and the pressure 
at step n + 1; 


n+l = jn+ l u * +u/j 

(3.109) 

t n-fl i ^ 

• = j p + pi ■ 

(3.110) 


The preprocessing stage of the solution needs to be done only once. The increase in 
solution time is very small in relation to the total time required to sufficiently advance the 
solutions to achieve a steady-periodic state. To solve 3.100, an Uzawa algorithm is used [46] 
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which decouples the pressure solve from the velocity solve and results in the two equations 


Sp* = DiH _1 B. (3.111) 

Hu* = dV + B1. (3.112) 

where S is defined as 

S =DiH- 1 D[ + D 2 H- 1 D^. (3.113) 

Equation 3.111 is solved with a nested conjugate gradient scheme and Equation 3.112 is 
solved for the velocity components with a conjugate gradient algorithm. 

3.2.2 System Solution Strategy — Iterative and Direct Solvers 

In this section, the solvers used to make the necessary matrix inversions at each step of the 
solution process are described. Matrices must be inverted for the solution of the preliminary 
Stokes solve (nested inversions for S in Equation 3.111 and of H in Equation 3.1 12). to obtain 
the approximate velocity solution in Equation 3.90. for the solution of the pressure correction 
term in Equation 3.91, and finally for the temperature solution update in Equation 3.96. 


3.2.2. 1 Iterative Solution — Conjugate Gradients 

The iterative solution algorithm of choice for finite element solvers is the conjugate gradient 
method. The conjugate gradient algorithm is extremely efficient for solving well conditioned 
problems and requires very low storage. The conjugate gradient algorithm has been exten- 
sively reported [31. 68] and the underlying theory is not detailed here. The preconditioned 
conjugate gradient algorithm for the solution of Ax = /, for a symmetric, positive-definite 
matrix A € IR nxn , preconditioning matrix M € IR nxn , x € JR n . and f € lR n is 

m = 0. r° = / - Ax°, 2 ° = w° = z°, 

while £ (r m ) < tol 
m = m 4- 1 

a m = {z m ,z m )/{w m , Aw m ) 
x m+1 = x m + a m w m 

r m + l _ r m _ Q m Au> m 
z m+l _ 

= (z m+1 .r m+1 )/{z n V m ) 

w m+l _ r m+lpm w m 

end 
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The stopping criterion is based on £(r m ) ? where 


e(r m ) 


n T m 


1 


Ik Bij 


(3.114) 


It should be noted that, for the preconditioning to be efficient, the solution 2 = M -1 r must 
be easily obtained. The preconditioner used for all of the iterative solves was simply the 
matrix diagonal M = [/l u ], i = 1, . . . , n which has only limited beneficial impact in practice. 


3. 2. 2. 2 Direct Solver — Sparse LU-Decomposition 

A sparse matrix solution package written by Kundert [43] has been used to perform the 
direct matrix solves when the Navier-Stokes code is run in one of the the direct-solver 
modes. The package is Sparse Version 1.3 and is available at Netlib. 1 The Sparse package 
is written in C and manages all of the memory allocation internally making it very easy to 
integrate into other codes. It can handle arbitrary square, real and complex, linear systems 
and is also able to find determinants and estimate ill-conditioning errors. According to its 
documentation, it performs as fast or faster than other sparse matrix packages and. at least 
anecdotally, it seems to be very efficient. 

The sparse solve consists of two stages. For the first stage, matrix factorization, the 
matrix elements are re-ordered and a LU factorization is performed [11], The lower and 
upper triangular matrices are then stored. The second stage, forward and backward sub- 
stitution, is performed on any number of right hand sides without having to re-factor the 
original matrix. For details of the algorithm, the reader is referred to [43]. 

As is well known, the factorization time is very large compared to each forward/backward 
solve. For the initial CPU cost of the LU factorization to be worthwhile, it must be prorated 
over enough time-steps to make the total solution time for the direct solve approach to 
be lower than for the iterative solver. The break-even point is easily exceeded by the 
total number of time-steps that are needed to obtain periodically varying, steady-periodic 
solutions. 


3.3 Sample Code Performance 

The Navier-Stokes code can be run in several modes, based on whether the iterative or 
direct solver is used for each matrix inversion. The lowest memory, most time consuming 
solution mode uses the iterative solver for all of the matrix inversions. Conversely, the 


'The Xetlib master index web address is http://netlib.bell-labs.com/netlib/master/index.html 
and the Sparse Version 1.3 package is located in the sparse directory 
http://netlib.bell-labs.com/netlib/sparse/index.html. 
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highest memory, fastest solution mode uses the direct solver for all of the matrix solutions. 
Various combinations of iterative/direct solves result in solution modes that fall in between 
the two extremes. The solution modes are described below. 

Mode I Conjugate gradient solvers used for all matrix inversions: H in Equation 

3.90, E in Equation 3.91. and H# in Equation 3.96. This mode requires the 
least memory but highest CPU time. 

Mode II Conjugate gradient solver used to invert H in Equation 3.90 and H# in 
Equation 3.96. The direct solver is used to invert E in Equation 3.91. 

Mode III Conjugate gradient solver used only to H# in Equation 3.96. The direct 
solver is used to invert E in Equation 3.91 and H in Equation 3.90. 

Mode IV Direct solver used for all of the required matrix inversions: H in Equation 

3.90. E in Equation 3.91. and in Equation 3.96. This mode has the highest 
memory requirement but the lowest CPU time. 

It is important to note that the H and H# matrices are dependent on the timestep size 
A t n . If the timestep size is changed during the time advancement, the H and H# matrices 
must be updated and. if running in Mode III or Mode IV. one or both factorizations 
must be recomputed. The time penalty incurred by frequent timestep adjustments and re- 
factorizations would rapidly overwhelm the performance benefit realized by the direct solver. 
To avoid such a situation, a conservative timestep is selected as soon as the periodic -steady 
state is achieved and is held fixed for the remainder of the time advancement. Because the 
E matrix is not A t n dependent, timestep updates can be performed more frequently when 
the code is used in Mode I or Mode II. 

To demonstrate the performance of the code, timings have been performed with the 
code running in each of the solution modes listed above. For the computational times that 
are reported, the factorization time has been neglected. It is reasonable to neglect this as 
the contribution to the total solution from the factorization for the 0(10.000) time steps 
that a solution is typically advanced is negligible. In all cases, the timestep size is constant 
throughout the run. 

The sample problem consists of a two-body, eddy-promoter configuration with the grid 
shown in Figure 3-2. The number of nodes and elements are indicated in the figure resulting 
in matrices sized as H € R 7 156x7 i 56 , E € jR 1842 * 1842 . and H# € JR 7156 * 7156 . Because the 
E operator is much smaller than H and Ho- and its inversion requires the majority of the 
fully iterative solution time, it is the best candidate for solution by the direct solver. This 
result is evident from the solution times and memory requirements listed in Table 3.1. 
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Figure 3-2: Mesh plot of the configuration used for code performance tests (N el = 3471, 
Af v = 7156. AfP = 1842). 


Solution Mode 

CPU seconds/timestep 

Core memory (Mbytes) 

Mode I 

7.91 

18.2 

Mode II 

2.51 

26.5 

Mode III 

1.32 

44.1 

Mode IV 

1.22 

63.4 


Table 3.1: CPU time and memory required for the four solution modes running on a Hewlett 
Packard C-160 workstation. 
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Chapter 4 


Baseline Surrogate Approach 


To begin the description of the baseline, non-parametric validated, surrogate framework it 
is first necessary to define a hierarchy of approximation models. Implicit in this discussion 
is the assumption that the governing equations (or experimental setup) accurately represent 
the problem under consideration. With this assumption, a hierarchy of three approximation 
models are defined: 

Exact This is the exact, analytical solution to the governing equations. For the eddy- 
promoter problem, it is the solution to Equations 2.25 — 2.36. The exact solution can 
be obtained only for the very simplified cases (e.g. channel flow with no obstructions). 

Truth The truth model is the 2 ncl -order in time. IP? — Pi finite element solution of the 
governing equations that is described in Chapter 3. While this is clearly not of the ac- 
curacy of the exact solution, it represents the best obtainable solution for the problem. 
The truth model serves as the reference for the validation and error analyses. 

Surrogate The surrogate model is the very inexpensive (by design) approximation to the 
exact solution and typically consists of a simple input-output model (e.g. response 
surfaces, scattered data interpolation, etc.). The surrogate model is used in place of 
the truth model in design studies. 

The non-parametric validated, baseline surrogate framework has been extensively re- 
ported [55. 57. 77, 78, 79] and is only repeated here to serve as a reference to the surrogate- 
Pareto optimization pursued in Chapter 5. In the baseline surrogate approach, a simple, 
inexpensive model (a surrogate) is substituted for an expensive simulation (or experiment as 
in Chapter 6) in the optimization process. The simulation serves offline first, to construct, 
and second, to validate the surrogates which are usually an inexpensive input-output func- 
tions. This is in contrast to online optimization approaches in which the simulation is called 
directly by the parent optimization process. In online approaches, the simulation serves as 
a subroutine call to the optimization routine. 
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Approximate, offline optimization approaches have been widely used but typically do 
not provide rigorous error estimates for surrogate-predicted designs. The distinguishing fea- 
ture of the non-parametric, validated surrogate approach is the surrogate validation step 
of the framework that provides a posteriori error estimates on designs near the surrogate- 
predicted. optimal design(s). The surrogate approach is related to probably-approximately- 
correct approaches [74, 26] and information-based complexity theory [73]. It differs, how- 
ever, from the former in that it is truly non-parametric (no assumption is made in regard to 
the distribution of optimal input points) and from the latter in that it requires no regularity 
estimates for the input-output function. 

The surrogate approach is best suited to problems that meet at least some of the follow- 
ing classifications: (1) The underlying simulation or analysis is computationally expensive 
making the large number simulation queries required for the optimization intractable. (2) 
The integration of the analysis code with the optimization package is impractical because 
of difficult interfacing or the lack of sensitivity derivatives. (3) The problem is global in 
nature and the design space covers a wide range of potential designs in which there may be 
many, locally-optimal designs. (4) Flexibility is needed in terms of rapid turnaround, as in 
cases in which the design goals evolve with the design and more is learned. 

In each section of this chapter, the surrogate framework is presented for a general prob- 
lem with two outputs and two performance metrics that are functions of each output as 
well of the input vector. An important distinction to be made between the general prob- 
lem presented here and those cases analyzed previously [55. 57. 77, 78, 79] is that here, 
surrogates for the simulation outputs are validated. The performance metrics which serve 
to quantify designer preferences are functions of the validated outputs as well as (possibly) 
of the inputs. This adds additional flexibility and generality to the surrogate framework in 
that the validated simulation outputs can be used in future design studies with different 
performance metrics and similar error bounds will obtain without the need for additional 
appeals to the simulation. In the earlier surrogate work, the performance metrics are val- 
idated directly. This restricts the design studies to that particular choice of performance 
metrics. A change in the form of the performance metrics therefore requires that the vali- 
dation step be repeated, necessitating additional appeals to the simulation. The two-body, 
eddy-promoter problem serves as an example of the two output, tw r o performance metric 
design problem presented in this chapter. 

The baseline surrogate framework consists of four steps. The first step of the framework, 
surrogate construction, is discussed in Section 4.1. In the second step of the framework, 
presented in Section 4.2, the inexpensive surrogates are validated against the truth simula- 
tion at randomly selected points. In Section 4.3 the third step of the framework, surrogate 
based design, is presented. The fourth and final step is the a posteriori error analysis of 
the surrogate-predicted designs and this is presented in Section 4.4. The a posteriori error 
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analysis, which uses the validation results without requiring additional appeals to the sim- 
ulation. gives error estimates for designs near the surrogate-predicted optimal designs as 
well as an estimate to the design optimality. 

4.1 Surrogate Construction 

The first step of the baseline surrogate framework is surrogate construction. To begin, given 
a vector of M m modeled inputs. p m G C R A ' m that lie in the model design (or input) 
space Q m . surrogates for two (K = 2) truth outputs of the simulation. <p[ (p^) : Cl m -► R, 
and (f>o{ p' m ) : — ► R- are constructed. Although the intent is that the surrogates, 

<MPrn) : -> R and (f>o ( p^ ) : R • approximate the simulation as closely as 

possible over <p\(p'm) ~ <t>\{v'm)' ^(Pm) ~ the results that will be obtained 

in the later sections of this chapter are valid regardless of the quality of the surrogates. 
The surrogates can be constructed by any means, whether it be by appeals to the truth or 
to another simulation, by empirical relationships, by the use of limiting solutions, or some 
combination of all of these. 

For the two-body, eddy-promoter problem, the two truth value outputs obtained from 
simulation are is the pressure forcing term required to achieve the target flow rate. i/-’o (p m ). 
and the inverse Nusselt number, #^(p m ). The surrogate construction set 

X‘° = { ( Pm • 9l (pm) - '0o" (Pm ) ) 1 (Pm-^O (Pm)- 4>o (P»> )) ,V co }• ("1 • 1 ) 

is formed by appealing to the finite element simulation (Truth) at 91 points in the design 
space fi,„ (defined in Section 2.6.1) and by including an additional 165 limiting solutions for 
a total constriction sample size of N co = 256 input-output pairs. The 91 simulation points 
have been chosen from an orthogonal array (Appendix D.2) that ensures a good space filling 
design. The additional 165 limiting solutions correspond to the upper and lower limits of 
Pm2 ( = 0 and p, n2 = 1) and represent duplicate geometric realizations at different 
points of the design space Q m . 

The approximation scheme used for the surrogates is a radial basis function fit through 
the construction set, X co , for each output. The particular radial basis function used here 
is described in Appendix D.l. The surrogates for the outputs are interpolating functions 
in that the output surface passes exactly through all of the points in the construction set 
X co . As demonstrated empirically in Section 2.5.4. the output functions are continuous 
over the entire model design space Q m and an interpolatory model is appropriate without 
the need to track and account for discontinuities. The plot in Figure 2-6 shows that the 
output function continuity is preserved across a topology change. Three-dimensional mesh 
plots of the surrogates are given in Figure 4-1. On the left in the figure are p ml -pm2 slices 
of the surrogates for p m3 = 0.50 and on the right are p ml -pm 3 slices of the surrogates for 
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Figure 4-1: Mesh plots ofthe output surrogates: (a) 0o(Pm)^slice for p m:i = 0.50. (b) 0o(P7n) 
slice for p m2 = 0.50. (c) ^o(Pm) slice for p m3 = 0.50. (d) ^o(Pm) slice for p m2 — 0.50. 


p m2 = 0.50. The temperature output 0o(Pm) is at the top and the pressure output ^o(Pm) 
is at the bottom. 


4.2 Surrogate Validation 

The second step of the surrogate framework is surrogate validation. In this step, the surro- 
gate models are validated against the truth at input points that are randomly chosen from 
the model design space f2 m . To proceed with the description of the validation step, several 
functions must first be introduced. 

The importance function p( p rn ) serves as a probability density function for the selection 
of validation points from the model design space Q m : 

I p(Pm) dp m = I- (4.2) 

■'Pm 

The importance function leads to the notion of a p-measure of a subdomain of the design 
space Q m which is simply the weighted relative A/ m -volume of the subdomain. For any 
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subdomain V C fl m . "the /^measure of V.'' p p {T>), is 

Hp{V) = [ p(p m ) dp m < 1. (4.3) 

Jv 

The p-measure of of full model input domain Cl m is one. p p (f2, n ) = 1. 

Next, K scaling functions are introduced, one corresponding to each of the output 
surrogates that are to be validated. For the general two output problem examined here, the 
scaling functions p 0l (p' m ) : fim dR- and Q^ip'm) • ~ ■+ ^R+ are the two strictly positive, 

error-scaling functions. The scaling functions serve primarily to make the errors between 
the truth value and the surrogate values for each of the outputs of the same magnitude. A 
secondary purpose for the scaling function is to reduce the impact of regions of the design 
space where large surrogate errors may exist. 

With the importance function and error-scaling functions defined, the validation sample 
set X va is formed as 


* la = {(P ml .^(Pml)^ r (Pml)).-.(Pm,V^[(Pm.v).^(P ra .v))}, P nl! ~ p(Pm)- 

(4.4) 

where the input points P r7U are independent, identically distributed random points drawn 
according the probability density function p(p m ). The truth outputs in 4.4 are obtained by 
appeals to the simulation at each validation point P mj . The size of the validation sample 
set is 


where \z] is the smelliest integer that is greater than 2 . Later, it will be shown that 
Si represents the p- measure of the uncharacterized region (the region in Q m where the 
surrogate error magnitudes |<?>{'(p' m ) — 0i(p' m )| and \(po (Pm) ~ are unknown) and 

e 2 is the significance of the probabilistic error estimates that are developed. 

With the validation sample set formed, the model prediction error U is computed as 


U = max max 

V ia \ 


'|tf(Pmi)-0l(Pmi)| |4(Pm.)-<MPr 


9o 1 (Pmi) 


9<?2(P'ru) 


(4.6) 


and the uncharacterized region T is defined as 


T = p m I max 


\4>T(P m) - jMjOl l^PPm) -02(Pm)l \ > y 


9q\ (Pm) 


9oi ( Pm ) 


(4.7) 


The model prediction error is the maximum, absolute truth-surrogate difference over the 

set of randomly drawn, validation input points P m ,.z = 1 N. The uncharacterized 

region is that subset of the model design space for which the scaled truth-surrogate 
error is greater than the model prediction error U . The validation provides no information 
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as to the location of the uncharacterized region nor does it provide any insight as to the 
magnitude of the errors in that region. This uncertainty manifests itself in the locality of 
the error estimates that are developed in Section 4.4. 

The model prediction error provides a probabilistic estimate as to the maximum truth- 
surrogate error over a fraction of the model design space The precise validation state- 
ment reads: 

Pr(/x p (T) < ei) > 1 — e 2 - (4.8) 

What 4.8 says in words is that, with probability of at least 1 — the fraction of the model 
design space for which the surrogate differs from the truth by more than U is less than 
€\. The relationship between N, , and £2 is entirely determined by the sampling theorem 
in Equation 4.5. The validity of 4.8 given 4.5 can be proved with order statistics. The proof 
is given in Appendix B.l. 

The validation step is the final step in the surrogate design framework where truth cal- 
culations are needed. The results of t lie validation framework provide precise, probabilistic 
bounds on the prediction error of the surrogates for designs near surrogate-selected opti- 
mal designs, without having to appeal to the simulation. The estimates can be generated 
for any number of designs and can be obtained very rapidly as the only required function 
evaluations are of the surrogates. This gives near-instantaneous turnaround for the er- 
ror estimates of each design and makes that design process interactive, which is extremely 
important in situations where the design goals evolve as more is learned. 

The validation step has been performed for the global surrogates of the two outputs for 
the eddy-promoter problem. Because no prior importance had been assigned to any region 
of the model design space Q m , a uniform probability density function p{p m ) = 1 has been 
chosen for the selection of the validation input points. A total of 24 simulation evaluations 
have been budgeted for the global validation which, from Equation 4.5, gives e 1 = 0.0561 
and £2 = 0.2500. The scaling functions have been set to constant values g$( p m ) = 2.00 and 
iMPm) = TOO. Forming the validation sample, X VQ . and finding the maximum, unsealed 
surrogate error for each output gives 

max|0(pp m ) - #o(Pm)| = 0.0968, (4.9) 

max|^(p m } -^o(Pm)| = 0.1910. (4.10) 

The model prediction error U is computed from Equation 4.6 and is U — 0.9681. From 
this, the following statement on the global quality of the surrogates can be made: With 
probability of at least 75% (1 - £ 2 ), 

\ej (Pm) - 0~o(Pm)| < 0.0968(7* (p m ) = 0.1936, (4.11) 
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(4.12) 


h/>o (Pm) - ^o(Pm)| < 0.0968g v (p m ) = 0.0968. 

over at least 94.39% (1 — eq) of the model design space. The global bound for the two 
output surrogates hold for the same 94.39%-sized subset of the model design space flm. 


4.3 Surrogate Based Design 


It is assumed that the design preferences are characterized by two performance metrics, each 
of which is a function of both surrogate outputs, <Mp' m ) : fi m ' -> iRand<fo(Pm) : ^m' Sti 
as well as of the full design vector p'. For both metrics, lower values are assumed to 
correspond to better designs. 

To form the surrogate metrics, the output surrogates. ~ </>f(Pm) and <fo( p' m ) ~ 

<p 2 (p( r , ) . are substituted into the performance metric functions. The resulting surrogate 
metrics $j = < Fl(</ ) l(p , m ); <fo(Pm)’ p ; ) : y ^ and ^2 = ^2(^l(Pm)- ^(Pm)- P*) : —■ y Hi 

are then used for the design studies. Because the surrogate outputs are trivial to evaluate, 
any number of designs may be pursued and the trade-offs between improvements in one 
metric versus the other can be fully examined and understood. The a posteriori error 
bounds developed in the following sections will apply for any design chosen based on the 
surrogate performance metrics and will require no additional simulation evaluations. In 
addition, because the surrogate outputs have been validated (and not the full performance 
metric functions), different performance metrics (using the same inputs, outputs and design 
spaces) can be used as the basis for design studies and the error estimates will still be 
derivable. This gives additional flexibility to the design process. 

The eddy-promoter problem that serves as the illustrative example here falls into the 
problem class just described. The output surrogates, constructed in Section 4.1 and val- 
idated in Section 4.2. are used to form the surrogate performance metrics. The output 
surrogates, 0o(Pm) and V’o(Pm)- constructed in Section 4.1 are substituted into Equations 
2.97 and 2.98. The resulting surrogate performance metrics are given by 


0(p) = Q{0 0 {pm)-.P) = log io 


1 

RePr 


+ #o(Pm) 


t)l\ ' 


(4.13) 


'I'(p) = ^(V’o(Pm).p) = log jo [v>o(p m )f?eV£] , (4.14) 

and are dependent on the full design vector p = (Pm^Pa) G Q = [0, lj '^ -4 . The range for 
the analytic input, m G [13.332, 106.656], is as stated in Section 2.6.1. 

For the design study, two designs, Q = 2, are pursued. The first design problem (q = 1) 
is to find the configuration that gives the best heat transfer. The second problem (q = 2) 
imposes a constraint on the value for the pumping power fl/. Both problems, q G Q = {1. 2}. 
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can be formulated as 


p 9 = arg min 0(^o(Pm),p)- Vg € Q (4.15) 

{pen | >{'(vo(pr7,).p)<u < '} 

where p 7 is the surrogate-predicted design point that minimizes the surrogate-based, wall 
temperature performance metric subject to the constraint that ty q = ^(^o(Pm)^ P 9 ) < ^ ■ 
The temperature performance metric at the optimal point p q is Q q = 0(#o(p£J, p 9 ). 

To perform the optimizations, a subspace searching, simplex algorithm (SUBPLEX) 
written by T. Rowan [61] and available on Netlib, 1 has been used. The algorithm is a 
generalization of the the Nelder-Mead [51] simplex algorithm in which high-dimensional 
problems are decomposed into low-dimensional subspaces that can be efficiently searched 
with the simplex algorithm. The SUBPLEX algorithm is well suited for the unconstrained 
optimization of noisy objective functions. The problem studied here (4.15) is a constrained 
optimization problem. The optimization results presented in this section have been obtained 
using the scalarization technique presented in Section 2.4.3. To find a global minimum for 
each design problem. 200 random restarts selected from a uniform distribution over Q have 
been performed. The inexpensive nature of the surrogates makes such a '’brute-force*’ 
approach tractable. The optimizer proved to be very robust in practice, and this was the 
primary reason for choosing it. 

The two design problems have been carried out for the eddy-promoter heat exchanger 
example, and the results are summarized in Table 4.1. For the first design, the pump- 
ing power constraint ip Q only needs to be set sufficiently high to ensure that it remains 
inactive. For the first design 4 } = 1 x lCr and for the second design 2 p = 9.4801. 

The optimization process gives the optimal points p 1 = (0.1008,0.2875.1.0000,1.0000), 
p 2 = (0.8549.0.4729,0.6449.0.0000), and the corresponding surrogate-predicted, perfor- 
mance metric values 0 1 = —2.2626, $ 1 = 12.5137, and 0 2 = —1.7106. 4> 2 = 9.4801. The 
inverse channel heights for each design are = 106.656, and t] 2 l == 13.332. The physical 
geometry for a single periodicity cell of each design point is given in Figure 4-2. 


Design ( q ) 


P 9 

0 9 


i 

1 x 10 ;3 

(0.1008, 0.2875. 1.0000. 1.0000) 

-2.2626 

12.5137 

2 

9.4801 

(0.8549. 0.4729. 0.6449, 0.0000) 

-1.7106 

9.4801 


Table 4.1: Baseline surrogate-based optimization results. 


^he Xetlib master index web address is http://netlib.bell-labs.com/netlib/master/index.html 
and the SUBPLEX software is located in the opt directory. 
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(a) Design 1. 
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(b) Design 2. 

Figure 4-2: Eddy-promoter configuration for each of the surrogate-predicted optimal de- 
signs. 

4.4 A Posteriori Error Analysis 

The a posteriori error analysis step of the surrogate framework takes as inputs the surro- 
gates. the probability density and error-scaling functions, the model prediction error, and 
the surrogate-predicted design points. From these inputs, and without appealing to the 
simulation, probabilistic estimates of the difference between the surrogate and the truth at 
points in the design space near the surrogate-predicted designs are derived. Two formu- 
lations for the estimates are presented: proximal region, in which the estimates are valid 
for a region of non-zero measure, and proximal candidate, in which the estimates apply to 
a randomly selected design point near the surrogate-predicted optimal point. In addition, 
the a posteriori error analysis provides an assessment of the surrogate-predicted design 
optimality in terms of an estimate of how much effort is required to improve the design 
beyond a given amount. The proximal region and proximal candidate error analysis meth- 
ods are presented in the next two sections, followed by the optimality analysis in the final 
subsection. 

4.4.1 Surrogate Predictability — Proximal Region 

To begin the Proximal Region error analysis, a model prediction region V* x € D is de- 
fined as a subdomain in the full design space Cl of p-measure ci, and that contains the 
surrogate-predicted design point p*. Here, the p-measure of V‘ x is computed only from 
the components in the model design space fil m . Recall that the design input vector p can 
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be decomposed into a modeled input vector p m (the inputs over which the surrogates are 
constructed) and an analytic component p Q (the inputs for which the performance met- 
ric responses are known analytically). Likewise, the surrogate predicted optimizer can be 
decomposed into p* = (pj^,p*). The prediction neighborhood can then be defined as the 
tensor product x p* where V^ nei is the model prediction neighborhood with 

the p-measure defined in Equation 4.3, fj, p ('P^ nsi ) = e The notion of how to precisely 
define the prediction region in practice is made more clear with the eddy-promoter 
example later in this section and in the following section. In general however, very limited 
assumptions are made as to the construction of V* x and it is not required to be connected. 

With the prediction neighborhood defined, probabilistic bounds on the outputs for de- 
signs near the surrogate predicted optimizer are first derived. The output bounds then serve 
to construct probabilistic bounds for the performance metrics for designs near the surrogate 
predicted designs. 

The predictability statement for the outputs reads: With probability at least 1 — gq. 


there will exist a region of non-zero measure r' m C V* m , i in the neighborhood of 
that for all points p' m E T' m 

lot < 4>J(p'm) < 

Pm such 
(4.16) 

l c , < 4>l( p'J < U 02 - 


(4.17) 

where the truth-output bounds are 



l 0l = min [<Mp' m ) -%(pL)]' 

{p'€V mzx } 


(4.18) 

u 0 i = . max [0i(p' m ) + Ug Px (p' m )}. 


(4.19) 

l P2 = min [fcip'm) ~ Ugo-,{p'm)}- 

{p'er- m;i } 


(4.20) 

u <t >2 = , , m a x J'MPm) +Ug< >2 (p , m )}. 
{p } 


(4.21) 


The above result effectively bounds the truth values of the outputs near the surrogate 
predicted designs. 

To extend the bounds to the performance metrics, a second step is pursued using the 
output bounds presented above. The predictability statement for the performance metrics 
reads: With probability at least 1 —£* 2 , there will exist a region of non-zero measure T f C 

in the neighborhood of p* such that for all points p' 6 T f 

< *i(*r(p' m ),4(Pm).P') < U* t , (4.22) 

L* 2 < $2(^f(p'J,4(Pm),P') < £/**, (4.23) 
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where the performance metric bounds are 


L<$, x = min $i( 2 i, 22 -p')t 

{p'G'P^ , 226 ^ 2 } 

(4.24) 

{ 7 $, = max ^ > i(^i,22 1 P , )t 

{p'€T>; i ,; l 62i,-'2€22} 

(4.25) 

L*, = min $2(21. 22- P^ 

{p , €P; l ,zi€ZuZ 2 €Z 2 } 

(4.26) 

U<po = max $2(21, Z-2-. v']- 

{p'ev; zieZuz-ieZi] 

(4.27) 


and the output ranges. Z\ and Z-i, are determined by Equations 4.18 — 4.21 and are 

Z l = {z\l 0l <z<u Ol }. (4.28) 

Z -2 = {z I /p., < z < Up.,}. (4.29) 


To define the prediction neighborhood Vt x precisely, a distance metric A(p ml .p m2 ) 
that quantifies the distance between two points is introduced. The prediction neighborhood 
around the surrogate predicted optimal point p? n is then defined as that region TZ. of size 
£q. that minimizes 


r A = max A(p m ,p)„). 

PmC'V 


(4.30) 


For the eddy-promoter problem, the distance metric function has been defined based on 
the outputs and is 


AfPrnn Pmi) = max 


^0 ( Pm 1 ) ^(Pmj) 


m 1 1 


^ofPrn-') 


(4.31) 


where hg and h v are positive scalars. The definition for the distance function in 4.31 can 
be used to specify the neighborhood with minimum 9 0 sensitivity (for h v sufficiently large) 
or. likewise, the neighborhood with minimum sensitivity (for hg sufficiently large). By 
setting the levels of the hg and h v to match the relative variation of the surrogate over the 
region 7 Z, a balance between the 9 0 and xp 0 sensitivity can be obtained. 


The output bounds in Equations 4.18 — 4.21 have been evaluated for each design. The 
output bounds have then been used to evaluate the bound quantities given in Equations 
4.24 — 4.27. The resulting predictability statements for the truth values of the performance 
metrics are given below. For Design 1, the statement reads: With probability of at least 
75%, there will exist points in T 1 C V] x such that for p' 6 T l 

11.1809 < $(V^ (p)J. p') < 12.8194. (4.32) 
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0 


Figure 4-3: Output plane plot of the two design point performance pairs (•) and the asso- 
ciated predictability boxes. 


-2.3899 < 0(^(p' m ) ; p') < -2.1152. (4.33) 

For Design 2. the statement reads: With probability of at least 75%. there will exist points 
in T 2 C V; { such that for p' € T 2 

8.4716 < *(^(p' m ),p') < 9.8613. (4.34) 

-2.3372 <0(^(p' m ),p') < -1.4601. (4.35) 

The predictability statements are valid for each design for the corresponding prediction 
neighborhood. For a given design q, the region T q C V q x for which the points satisfy the 
performance metric bounds is identical for each metric. 

The results given above have been plotted in the output plane in Figure 4-3. In the 
figure, the surrogate performance metric values for each design are plotted as solid dots and 
the bounds on the truth performance metrics are plotted as boxes. It is evident from the 
predictability bounds for the two eddy-promoter designs given in 4.32 — 4.35 that the global 
output surrogates are quite poor. Even for the two very different designs that have been 
selected (based both on the associated geometries and the surrogate-predicted performance 
of each), the error bounds between the two design overlap. This implies that the surrogates 
are not accurate enough to distinguish between the two designs as regards temperature. 
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4.4.2 Surrogate Predictability — Proximal Candidate 


The Proximal Candidate a posteriori error analysis requires that we define a prediction 
region V * C Q that is of p-measure cr, and that contains the surrogate-predicted optimal 
point p*, in a manner similar to that presented for the Proximal Region analysis of Section 
4.4.1. A sample candidate design P^, is selected randomly according to ~ pp{p m ). The 
density function pr(p m ) is defined as 


PviPm) — P(Pm) 

a 


Vp m € V„ 


(4.36) 


where V^ na C is the sub-manifold of V* that corresponds to the model design space 
and the full neighborhood is the tensor product V x p*}. The full candidate 
design vector P* is given as P* = (P^.p*) where p* is the analytic inputs of the sur- 
rogate predicted optimizer p* = (p^.p*). The candidate designs can easily be selected 
with acceptance-rejection techniques without having to formally construct the prediction 
neighborhood V* ng . 

A second small parameter, c c . is introduced and is related to p-measure of V' na , cr. by 


1 


a(N+ 1) 


(1 - (1 -cr) 


•V— 1 1 


(4.37) 


where 0 < s c ,cr < 1. and N is the validation sample size given in Equation 4.5. It will 
become evident below that s c is the significance of the probabilistic bound estimates that 


are derived. 


Now, given the two inputs to the analysis e c and cr. and the model prediction error U. 
the following statement can be made: With probability of at least 1 — e c . the truth output 
values at PJ), are bounded by 


where the bound values are 



(4.38) 

1%, < 4( Pm) < 

(4.39) 

Ipn 

II 

J! 

►d> 

i 

c: 

o 

3J 

(4.40) 

= $i(P* m ) + Ug 0i (P* m ), 

(4.41) 

ll 2 =fo(P; n )-Ug 0 . 2 (P* m ). 

(4.42) 

4 = ^(p;) + %(p;). 

(4.43) 


Equations 4.38 — 4.43 bound the truth value of the outputs for the candidate design point 
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To extend the analysis to the performance metrics, a procedure similar to that of Section 
4.4.1 is pursued. Again, the following statement can be made: With probability of at least 
1 — £ c , the truth performance metric values at are bounded by 



(4.44) 

l% 2 < $ 2 (^(P^),4(P^),P*) < u% 2 . 

(4.45) 

where the performance metric bounds are 


L% = min $i(*i, *2, P*)- 

1 {z^Z\,z 2 ^Z c 2 } 

(4.46) 

UL = max $i(2 1 .2 2 -P*)- 

1 {z,€Z\,Z2€Z%} 

(4.47) 

L% = min $ 2 (zi,.z 2 . P"), 

2 {*,€Zf,sa€ZS} 

(4.48) 

= max $ 2 (2 i,2 2 ,P*), 

{.*1 €2^:262!} 

(4.49) 

and the ranges for the outputs, Z\ and at the candidate design point 

Equations 4.40 — 4.43 and are 

are obtained from 

z\ = {z\r 0x <z< upj, 

(4.50) 

Z 2 C = 1 l c Q2 < z < U% 2 }. 

(4.51) 


The proximal candidate error analysis produces estimates that apply to a specific design; 
the randomly chosen candidate design P* near the surrogate-predicted optimal point p*. 
This is in contrast to the proximal region analysis presented in the previous section in which 
the error estimates applied to a region of finite measure in the neighborhood of p*. The 
proximal region analysis assures that the region will exist with 1 — €2 confidence, but can 
not provide precise information as to where in V* x the region exists. Although the proximal 
region analysis does provide a sense of stability, it may not be satisfying in some design 
scenarios. The error estimates for the candidate design developed in this section provide 
an alternative interpretation. In this scenario, the designer can be assured that the error 
bounds apply, with 1 — e c confidence, to the specific, randomly drawn candidate design. 

For the eddy-promoter problem, the same two surrogate-predicted designs given in 
Table 4.1 are considered. The uncertainty parameters, a and e c , in Equation 4.37 are set to 
a = 0.1578, and e c = 0.2500 which agree with the N = 24 global validation points that have 
been computed in Section 4.2. A new distance metric that quantifies the distance between 
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(4.52) 


two input points. p m i and p m2 , is introduced to define the prediction neighborhood 

A c (p ml , Pm 2 ) = [(Pml — Pra')) ' (Pml — Pm' 2 )] 2 ' 

where the [() ■ ()] term above is the vector dot product. The distance metric defined in 
Equation 4.52 is simply the Euclidean distance between the two input points in the model 
input space Q m . The metric in Equation 4.52 emphasizes designs that are "near" the 
surrogate— predicted optimal designs based on the eddy-promoter geometry (model inputs) 
whereas the distance metric given Equation 4.31 defines designs “near’ the surrogate- 
predicted optimal design as similarly performing designs. The prediction neighborhood is 
then defined as that region H, of size a. that minimizes 

r^ c = max A c (p m .p? n ). (4.53) 

pm€'/C 

With the prediction neighborhood defined as above. Monte-Carlo sampling is used to em- 
pirically measure the neighborhood and acceptance -reject ion techniques aie used to obtain 
a candidate for each design. 

For Design 1. the candidate design input vector is drawn according to the density in 
Equation 4.36 is P 1 = (0.225C. 0.4976. 0.821 1. 1.0000) and the corresponding surrogate per- 
formance metric values are 0* = 0(0 o (P, l J. P l ) = -2.2192 and = 

12.6431. The predictability statement for the candidate design point reads: With confidence 
of at least 75%. the truth performance metric bounds for the candidate design are bounded 

by 

-2.3745 < ©(^(P^J.P 1 ) < -2.1050. (4.54) 

12.4085 < < 12.7946. (4.55) 

For Design 2. the candidate design input vector is P 2 = (0.7928.0.3888.0.6837,0.0000) 
and the corresponding surrogate performance metric values are © 2 = ©(#o(Pm)-P 2 ) = 
-1.6766 and 'F 2 = < &{'ipo(P 2 ), Pm) — 9-4762. The predictability statement for the candidate 
point for Design 2 reads: With confidence of at least 75%, the truth performance metric 
bounds for the candidate design are bounded by 

-2.3372 < 0(^(P^),P 2 ) < -1.4601. (4.56) 

8.4716 < tf(V><T(P 2 ),Pm) < 9-8613. (4-57) 

The two design points and candidate surrogate performance metric pairs, and the pre- 
dictability boxes are plotted in the surrogate performance metric plane in Figure 4-4. The 
predictability boxes show graphically where the truth performance of the candidate design 
(shown as an open circle) will lie with the associated 75% confidence. Comparing the re- 
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0 

Figure 4-4: Output plane plot of the two design point peformance pairs (•), candidate 
design performance pairs (o), and the predictability boxes for the candidate designs. 


suits in Figure 4-4 with the proximal region analysis results plotted in Figure 4-3 show 
that the predictability boxes for the candidate designs are sharper than for the proximal 
region. However, although both analyses indicate that the surrogates can discriminate be- 
tween the designs based on the pumping power metric, neither can do so in terms of the 
temperature performance metric. This further reinforces the conclusion that the surrogates 
are not performing adequately to select designs based solely on the surrogate predictions. 
The geometry for each of the candidate designs is plotted in Figure 4-5 and these can be 
compared to the surrogate-predicted optimal designs shown in Figure 4-2. 

4.4.3 Design Optimality 

In earlier work [77, 78, 79], optimality estimates were obtained for surrogate-predicted 
optimizers in the neighborhood that relied on the assumption of quasi-convexity of 

the truth response in the region The resulting optimality bounds, however, are often 

times too large to be of use. Another approach to the optimality analysis is to estimate 
the effort required to improve upon the surrogate-predicted design. The procedure is to 
develop estimates as to the amount of improvement in each metric that could be realized by 
randomly drawing additional simulation points according to the probability density function 
used for the validation step. It is again assumed that lower values of the two performance 
metrics are preferred. 
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(a) Candidate design 1. 

loo i 

I I 

(b) Candidate design 2. 

Figure 4-5: Eddy-promoter candidate design configurations for each of the surrogate- 
predicted optimal designs. 


Given a surrogate predicted optimizer p* with associated surrogate performance metric 

values = $i(0i(p; n ).^2(Pm).P*) and $*, = $2(<MPm)- <?’(Pm)- P*)- the exteilt that 
a surrogate-predicted, optimal design, p*. can be expected to be improved by randomly 
selecting design points, and evaluating the truth performance metric values for each, is 
bounded. For the analysis of this section, it is assumed that, because lower values of the 
performance metrics correspond to preferred designs, a lower value in at least one metric 
without an increase in the second, represents an improved design. This assumption has been 
shown, in Section 2.4.2. to be consistent with a wide class of multiobjective optimization 
formulations. 

First, an unbounded set of points in the surrogate, performance-metric plane is defined 
as 


A + = {s e JR 2 |si < $l(<MPm)' < MPm)-P*)’ *2<^l(Pm)-^(Pm)-P*)} ( 4 - 58 ) 

which is the set of all performance metric pairs to the left, and below, the surrogate- 
predicted optimal point performance. As stated in the preceding paragraph, input points 
with performance metric pairs in can arguably be called better than p*. A second set 
of points in the surrogate, performance-metric space is defined as 

A 0 pt = {s € JR 2 | 3p' € O', z i € Z° x pt , z 2 € ZT s.t. 

^ 1 (01 (Pm) • ^(Pm)' PO < 5i.$2(0l(Pm)-02(Pm)-P') $ ( 4 - 59 ) 
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Figure 4-6: Surrogate performance metric plane plot showing the surrogate predicted opti- 
mizer performance (•). the expanded optimality set A opt (light shade), the uncharacterized 
optimality set ^4 unch (dark shade), and boundary dA opt . 


where the ranges and Z 9 Pt are 

Z° pt = {z I Vp' m G P.' m . [4>\{v',n) - ^Oi(p'm)] < X < [<MPm) + ^Pi(Pm)]}- (4.60) 

ZT = {2 I Vp' m G f>' m , [^(p' m ) - Ug 0i (p' m )] < z < [foip'm) + ^5p 3 (Pm)]}- ( 4 -61) 

The third region in the surrogate performance metric plane is then simply defined as *4 unch = 
A + \ *4 0pt . The region A unch is the set of possible performance metric values that improve 
upon the surrogate-predicted design performance and correspond to input points in the 
region of the design input space that is uncharacterized by the surrogate validation. A 
schematic of the region ^4 0pt (shaded) and >T inch in the surrogate performance metric plane 
is given in Figure 4-6. In the figure, the surrogate performance metric response pair of the 
surrogate predicted optimizer plotted as a solid dot. 

Next, a sequence of N opt , independent, identically distributed, randomly selected inputs. 
P mi . . . . , P m vopt , is drawn according to the probability density function p(p m ); P mt ~ 

p(p m ).i = 1 , N opt . The optimality statement then follows that, the probability that at 

least one of the randomly selected design points will have truth performance values that lie 
in .A unch is less than 1 — The number of random input vectors drawn, jV opt , is computed 
as 

N° p* < N (^- - l) , (4.62) 

where the parameter Si E]0, 1[ is set independently. The optimality result states that, with 
a confidence of less than 1 — ei. N opt additional truth evaluations at randomly selected 
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points will produce a design with truth performance metric values better than *4 opt . 

For the eddy-promoter design problem, the upper boundaries for Af have been 
plotted as dashed lines in Figure 4-7. The boundaries of ,4° pt and are plotted as thick 
solid lines and the uncharacterized optimality regions, >l“ nch and .4.2 nch , are shown shaded. 
Note that the regions ^i nch and ^2 nch overlap in the lower left corner of the plot even 
though the lightly shaded *42 inch supersedes ^4]* nch in that region. 

With the optimality regions evaluated, the optimality assessment follows directly. First, 
for Design 1, if is set to Si = 0.25, then according to 4.62, if an additional N opi = 72 
input points are drawn randomly, and the truth performance metrics for each evaluated, 
then the probability that at least one of the randomly chosen inputs, i, has truth metric 
output pair such that 

(Wo'(Pm,),Pi),0(fio'(Pm 1 ),Pi)) eAT C \ (4.63) 

is less than 1 -s L = 75%. Similarly, for Design 2. if an additional N opt = 72 input points are 
drawn randomly, and the truth performance metrics for each evaluated, then the probability 
that at least one of the randomly chosen inputs, i. has truth metric output pair such that 

(«(^(P m ,),Pi).©(^(Pm I -).P,-)) eAT ch , (4-64) 

is less than 1 — si = 75%. 

The optimality analysis provides an additional evaluation as to the quality of the sur- 
rogates. For very accurate surrogates, which would in turn result in an extremely small 
model prediction error U , the boundary on the regions >f^ nch and would be very close 

to the surrogate-predicted optimal performance pairs. This would suggest that even with a 
great deal of computational effort, the performance of a given design can be improved only 
slightly (and only at best with a 1 — £i certainty) and is likely not worthwhile. However, 
the model prediction error for the eddy-promoter problem is fairly large and this is evident 
by the large improvement that can be realized (the large distance between *4\ mch and ^4-2 nch 
and the design point performance pairs) for the additional computational effort. The ad- 
ditional effort, would likely be worthwhile for the eddy-promoter problem suggesting that 
improvements in the surrogates are required. 

4.5 Baseline Surrogate Summary 

The primary drawback to the baseline surrogate approach is the difficult construction and 
validation of the surrogate in high dimensional input spaces. This can be easily illustrated 
by considering the uniform importance function p(p) and a neighborhood of p-measure 
Si in the input domain Q = [0, l]* u . The neighborhood will span at least ;| , U in one of 
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Figure 4-7: Uncharacterized optimality regions shown shaded for each of the surrogate- 
predicted designs (•). 


the input directions which rapidly approaches one as M — > oc. The loss of localization 
as M -*> oc produces a corresponding loss in predictability through increased value of the 
neighborhood radius as given in Equation 4.30. 

In certain instances, the surrogate approach can be effectively applied to problems with 
high dimensional input spaces. One such case occurs when the shape inputs are highly cor- 
related. An example is the optimization of smooth body profiles in which highly oscillatory 
geometries are not likely optimizers [55, 56. 57]. For cases in which the inputs are highly 
correlated, the effective input dimension of the problem is reduced. The eddy-promoter 
problem is a second instance in which the surrogates can be effectively applied to a prob- 
lem with a high input dimension. By pursuing a Pareto formulation for the two-criteria 
optimization problem, the effective dimension for the validation and predictability error 
analysis is reduced from the number of inputs to one fewer than the number of critera: here 
a single parameter. The surrogate-Pareto formulation has been applied previously [36]. In 
Chapter 5, the surrogate-Pareto formulation is applied to the eddy-promoter problem. 
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Chapter 5 


Surrogate— Pareto Approach 


In this chapter, the nonparametric-validated, surrogate-Pareto optimization framework is 
presented. The basic structure of this chapter follows the presentation employed to describe 
the baseline surrogate framework in Chapter 4. The identical Exact/Truth/Surrogate hier- 
archy of solution approximations will be assumed. 

In Chapter 4, the weaknesses of the baseline surrogate framework in terms of design 
predictability was apparent for problems with more than "several" inputs. In fact, even 
for the eddy-promoter problem with only 3 modeled inputs, the surrogate construction is 
quite difficult due to the complexity of the input-output relationship. Although pursuing a 
surrogate-Pareto approach to the design optimization can not avoid the difficulties of sur- 
rogate construction, it does dramatically improve the predictability of designs for problems 
with many inputs and 2 or 3 performance metrics. 

In the following sections, the steps of the surrogate-Pareto framework are presented. The 
format for each chapter is to present the theory for a general, two output, two performance 
metric problem and then to illustrate the techniques with the eddy-promoter problem which 
is a specialized case of the general problem. Throughout the chapter, it is assumed that 
lower values of the performance metrics are preferred although no such assumption is made 
for the outputs. In the first section, the surrogate construction stage is presented. Second, 
the validation stage of the framework is presented. Third, the surrogate-based design step 
is described. Fourth, the a posteriori error analysis is described. 

5.1 Surrogate Construction 

Surrogate construction takes two stages. In the first stage, global surrogates are constructed 
to approximate the truth input-output function over the model design space Q m . The sur- 
rogate models are i/^(p m ) ~ ^o(Pm) and 0 o (Pm) ~ ^o(Pm) and are constructed exactly as 
they were for the baseline surrogate framework in Section 4.1. In the second stage of the 
surrogate construction, the global model surrogates are used to form the surrogate perfor- 
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mance metrics, which are then used to find the surrogate Pareto-optimal input manifold dA 
and output manifold C A . The construction stages are presented in the next two sections. 

5.1.1 Global Surrogates 

The global surrogate construction step is identical to that presented in Section 4.1. Given 
M m modeled inputs and the associated input (design) vector p m € C iR A/m that 
lies in the model design space Q m , two truth outputs ( K = 2), </>f(p' m ) : -> 1R and 

< ^2 (Pm) : ft m — > JR. are considered. Surrogates for the truth outputs are constructed, 
<P\{Pm) - 4>\ (pm)’ <fe(Pm) % 02 (P m)? ^ at approximate the truth input-output function 
as closely as possible. The surrogates can have any form, incorporate data from any source, 
and in general, no limiting assumptions are made concerning construction or quality. The 
validity of the results obtained later are not impacted by the surrogate approximation 
quality, however the utility of the results will be. 

The global eddy-promoter surrogate models. 0 o (Pm) and constructed in this 

step of the surrogate-Pareto framework are identical to those constructed for the baseline 
surrogate framework (Chapter 4). The surrogates are a radial basis function fit to N co = 256 
construction input-output pairs. The radial basis function is described in Appendix D.l. Of 
the 256 construction points. 91 have been obtained by appealing to the truth simulation at 
model input points selected based on an orthogonal array (Appendix D.2). The remaining 
165 construction points are duplicate points in the design space fi m that arise from the 
mapping from geometric inputs Z ep to the normalized input vector p m described in Section 
2.5.3. The duplicate points correspond to input points for which p m2 = 0. 1. Slices of the 
surrogate surfaces are plotted in Figure 4-1. 

5.1.2 Surrogate PO Manifolds 

The goal of the design problem is find the design (or designs) that achieves lower val- 
ues for the two performance metrics, (pJJ.p') and (Pm)>P')- 

Because the truth values of the outputs are typically very expensive to evaluate, the 
global surrogates, </>i(p' m ) and <fo(Pm)’ constructed in the previous section are used to 
form the surrogate performance metrics, $i(p) = ^i(^i(p , m ).^2(p , m ),p') : fi' -> R and 
3* 2 (p) = *$2(01 (Pm)’ foip'm)- pO : ft ; — » 2R- Because the output surrogates are (perforce) 
computationally inexpensive, the surrogate Pareto-optimal manifolds can be obtained to a 
sufficient resolution within a reasonable amount of time. 

To determine the surrogate Pareto-optimal input manifold C A , the scalarization proce- 
dure described in Section 2.4.3 is used. With this approach, the multicriteria optimization 
based problem on the performance metrics $i(p) and ^(p) is reduced to a series of scalar 
problems parameterized by the scalarization parameter w G W = [0, 1]. A min-max for- 
mulation [22] for the scalarization has been used as it ensures that the full C A is obtained. 
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The parameterized, scalar optimization problem has the form 


£(w) = arg min max(u;<I>t (p), (1 — u>)$ 2 (p)). w € VV/,, 


p'efJ' 


(5.1) 


and is solved for a sufficiently resolved, discrete set W/, = {uq, . . . , w$po } C VV. The 
curve £(u/), with duplicate points, and horizontal and vertical segments removed, is the 
surrogate Pareto-optimal input manifold £ A . For each value of w , the solution to Equation 
5.1 requires a global optimization. If the computational cost of or $2 is even moderately 
high, the total computational effort required to obtain C A would be prohibitive. 

With the surrogate Pareto-optimal input manifold. C A . known, the surrogate Pareto- 
optimal output manifold dA can be determined from the surrogate performance metric 
functions. The surrogate Pareto-optimal output manifold is 


M = ($i(((w))i,(((w))). 


(5.2) 


The result from the above process is similar to the curves shown in Figure 2-3 but. because 
the surrogates are only approximations of the truth outputs, there is no assurance that dA 
and C A will correspond to dA and C A . respectively. The validation and error analysis steps 
of the surrogate Pareto framework presented in latter sections of this chapter address the 
impact, that the surrogate errors have on the results. 

The surrogate performance metrics for the eddy-promoter problem are obtained by 
substituting the surrogate output functions constructed in Section 5.1.1 in place of the 
exact output functions in Equations 2.98 and 2.97. The surrogate metrics, 


0(p) = 0(<9o(Pm).p) = log 10 


1 

RePr 


+ #o(p m) — ^ 
VL ' 


(5.3) 


$ (p) = ^o(Pm).p) = lo §10 [t/>o(p m)-Re 3 pl] , ( 5 . 4 ) 

are dependent on the full design vector p = (Pm-Pa) G = [0. The range for the 
analytic input, tjl 6 [13.332.106.656]. is as stated in Section 2.6.1 and the value of rji is 
prescribed by the design input p a = p\ G [0, 1] and is 


HL dLmin P-i(VLma\ Pt-min) - (5-5) 

The next step of the process is to define the surrogate, performance-metric achievable set 
as 

A = {s G R 2 | 3p E 9. s.t. 0(p) < si, $(p) < S 2 } (5.6) 

and find dA. the boundary of A not at infinity, with the min-max scalarization procedure 
in Equation 5.1. The surrogate performance metrics are used and the sequence of scalar 
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optimization problems that must be solved is 

£(w) — argminmax(u>0(p). (1 - lu)'I'(p)). w G W/i (5.7) 

where C W = [0,1]. Because the surrogates are (by design) trivial to evaluate, solv- 
ing for £(W/i) is a reasonable computational task. Once ^(W/J is obtained through the 
scalarization process, the surrogate PO output manifold is 

dA = (©of (W*)), Mmh))). (5.8) 

Intervals of [0, 1] that are not PO are assumed removed and W is appropriately rescaled so 
that the full extent of W = [0, 1] is PO. 

To solve the optimization problem in Equation 5.7. a subspace searching, simplex opti- 
mization package, SUBPLEX, has been used [61]. The SUBPLEX method is based on the 
Nelder-Mead simplex algorithm [51], but reduces the searches to lower dimensional mani- 
folds to more efficiently solve problems with high input dimensions. The primary reasons 
for the using the SUBPLEX package so solve the scalar optimization problems its demon- 
strated robustness, its availability at Netlib, and its ease of integration with the radial basis 
function surrogate routines. Because the SUBPLEX algorithm converges upon local opti- 
mizers. random restarts were used to ensure that a global optimal point was obtained for 
each value of w. " 

The surrogate PO manifolds have been obtained by solving Equation 5.7 over a suffi- 
ciently refined, discrete set WV The surrogate PO output manifold is plotted as a solid line 
in Figure 5-1 and the corresponding PO input manifold is given in Figure 5-2 in four plots: 
in each plot, one component of the input vector p has been plotted versus the scalarization 
parameter w. An inner nested, bisection approach has been used to isolate points in VV^ for 
which the PO input manifold is discontinuous. In total, Equation 5.7 has been solved for 
N PO = 932 values of w to obtain sufficient data to reconstruct the PO manifolds accurately 
with linear interpolation between points in VYV In the presentation of the results, it is 
assumed that V\4 ss W to sufficient accuracy to ignore the discrepancies. Even with the 
very inexpensive radial basis function surrogates, the solution to the scalarization problem 
at 200 u;-points, with 300 random restarts for each value of w. requires approximately 3.6 
CPU hours on a Hewlett-Packard C-160 workstation. 

The surrogate PO output manifold. <$L4. for the eddy-promoter problem is plotted 
in Figure 5-1 as a solid line. The corresponding curve for a plane Poiseuille, heat ex- 
changer (Appendix A) is labeled dA pp and plotted in the figure as a dashed line. The 
only input available for the plane-channel heat exchanger is rj pp and the range required 
to achieve the same range of temperature performance metric as the eddy-promoter is 
, pp = [24.546,507.849]. The third line (plotted with a dash-dot) in the Figure is the PO 


98 




Figure 5-1: Surrogate-predicted. Pareto-optimal output manifolds for the eddy-promoter 
heat exchanger. OA. the plane Poiseuille heat exchanger. 0A PP . and for the construction 
points. dA co . 


output manifold for the N co = 256 construction points. dA ro . It is obtained by evaluating 
the performance metrics for each construction- point, truth-output value over the range 
Hl G [13.332. 106.656]. forming the construction achievable set. and then finding the points 
that are PO. 

The PO output manifolds plotted in Figure 5-1 suggest that a considerable performance 
improvement is realized for the eddy-promoter exchanger relative to the simple channel- 
flow exchanger. This is evident by considering the required pumping power. 'P. for each 
exchanger for a given temperature performance value. 0. In all cases, to achieve a given heat 
transfer performance, the plane-channel exchanger requires a higher pumping power, and 
hence, is less efficient. The comparison between dA and dA pp shows that the surrogates at 
least predict better performance than obtained by simply selecting one of the construction 
point configurations. The uncertainty in the accuracy of the output surrogates, however, 
calls into question whether the eddy-promoter heat exchanger would actually perform better 
than the plane-channel exchanger, and the full set of construction point configurations, in 
practice. The surrogate validation and a posteriori error analysis steps will address this 
question. 
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Figure 5-2: Surrogate-predicted. Pareto-optimal input manifold plotted for each input 
p E ft versus the scalarization parameter w . 


5.2 Surrogate Validation 

The validation step of the surrogate-Pareto framework requires two steps. First, the sur- 
rogates are validated over the full input domain ft as was done in the validation step of the 
baseline surrogate framework. The results of the global validation are used in the optimal- 
ity analysis. Second, the surrogates are validated in the lower dimensional manifold of the 
design space "near* the PO input manifold. The results of the second validation step are 
used for predictability analyses of surrogate selected designs. 


5.2.1 Global Validation 

The global validation step for the surrogate Pareto framework is identical to the validation 
step of the baseline surrogate framework. The output surrogates 0i(p' m ) and faiPm) are 
validated over the model input domain ft' m . The details of the global validation step are 
given in Section 4.2 and are only briefly repeated here. 

To begin, a sequence of N validation points. P mi , . . . , P m . v , is randomly drawn accord- 
ing to the probability density function p( p m ) : ft' m — > JR, P mi ~ p(p m ). The truth outputs 
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are evaluated at each validation input point and the validation sample set X va is formed 


X va 


{(Pm 1 ^[(Pml).^(Pml)),...(P« i V^[(P«A-)^(Pm i v))} 


(5.9) 


The size of the validation sample set is exactly that given in Equation 4.5 and is 


N = 


In £2 

ln(l - e t ) 


(5.10) 


where [ 2 ] is the smallest integer that is greater than 2 . 

The global model prediction error U is computed by finding the maximum scaled differ- 
ence between the surrogate and the truth over the full set of validation input points. The 
model prediction error is 


U = max 

x va 


l<tf(P r 


max 


0l(Prm)l |02*(P mi) 4>2 (Pr 


9o\ (Pr 


<?£> (Pmi) 


(5.11) 


where g 0l ( p' m ) : Q' m -4 R~ and g 0 ,( p' m ) : -> R- are the two strictly positive, scaling 

functions. Introducing the notion of the uncharacterized region T as defined in Equation 
4.7, the following validation error statement can be made 


Pr(//. p ( T) < £ 1 ) > 1 - f 2 - (5.12) 

The statement in 5.12 gives an assessment of the surrogate performance over the model 
input domain fi' m . It bounds, with a confidence of 1 — £ 2 , the maximum scaled difference 
between the surrogates and the truth over 1 -5i of the design space. The proof of the result 
given in Equation 5.12 given the validation sampling theorem of 5.10 is given in Appendix 
B.l. 

The global validation step for the eddy-promoter problem follows exactly the validation 
step in Section 4.2. The validation sample size has been set to N = 24 and the validation 
input points have been selected according to a uniform probability density function p(p m ) = 
1. From the sampling theorem in Equation 5.10, E\ — 0.0561 and e -2 — 0.2500. The scaling 
functions have been set to constant values ge{p m ) = 2.00 and 5v(Pm) = 100. The validation 
sample set has been formed and the resulting model prediction error is U = 0.0968. This 
result is then used in Section 5.4.3 for the optimality analysis. 


5.2.2 Surrogate PO Design Validation 

The goal of the surrogate PO validation step is to provide information that will be used for 
predictability analyses of surrogate-predicted optimal designs. The strategy is similar to the 
global validation described in section 5.2.1. but instead of validating over the entire model 
design space Q m . the validation points are restricted to the vicinity of the surrogate PO 
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A 


Pm l 

Figure 5-3: A schematic the Pareto-optimal input manifold (solid line), the PO design 
space (shaded), and the model design space fi m ' (dashed lines). 



input manifold C A . As will be seen later, this restriction greatly improves the predictability 
of the a posteriori error analysis results. The improved predictability translates directly 
to error statements that are much more sharp, and hence, more meaningful and useful. 
As was demonstrated in section 2.4.2. the strategy outlined here achieves the improved 
predictability without losing design generality and flexibility. 

For the general problem of validating the surrogate outputs. 0i(p' m ) : Q f m — ► JR and 
02 (Pm) : Q'm JR- in the vicinity of the PO input manifold C A , a PO design space 
c fi' m is introduced. The PO design space is a finite-width (but narrow) i4 tube“ 
in the model design space S2' m that contains the surrogate PO input manifold C A . A 
schematic of is shown shaded in Figure 5-3 for a two-dimensional model design space 
Q' m . In the figure, the PO input manifold is plotted as a solid line and the boundaries of 
n' m is shown as a dashed line. 

An importance function, analogous to the importance function p( p m ) used in Section 
5.2.1, is introduced pc{ p ; m ) * JR- The PO importance function serves as a probability 

density function for the selection of validation points in It is also used to introduce 

the notion of the p-measure of a set 71 C 

Pp(K) = [ PciPm) d P m < L (5-13) 

J7Z 

which is the weighted volume of 7Z relative to the PO design space note p p (VL^) — 1. 

Because of the complexity of the function that describes C A , it is often times not possible, 
or practical, to attempt to define a particular density function. pc{Pm)i explicitly. Instead. 
Pc{p'm) defined implicitly through the procedure used to select random, validation input 
points P ~ PciPm)- This approach is straightforward to describe and easy to implement in 
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practice. First, a random scalarization parameter is introduced 

W ~ fw{ w ') (5-14) 

distributed according to the probability density function f\v{w'). The scalarization param- 
eter specifies the coordinate of a point along the PO input and output manifolds. It can be, 
but is not restricted to, the scalarization parameter w used to obtain the PO manifolds in 
Section 5.1.2. Another coordinate such as the one of the performance metric values along 
the PO output manifold [36] may be used as well. Next, a uniformly distributed, unit vector 
V' and a radius ry> are introduced. The density function is that function, Pc(p' m )i t ^ iat 
prescribes the distribution of the validation points 


P'rn ~ Pc(p'm)- 


(5.15) 


such that the random input points are given by 

P',n =Z(W , ) + r v ,V'. 


(5.16) 


and the validation sample set is 


- { Pm e 


Ip m 


— <f(u/)]| < r\ >, Vu/ € w'}. 


(5.17) 


Although this implicit definition for pc(Pm) is difficult to express analytically, it is easy to 
sample in practice. 

With the probability density function pc{p' m ) defined as above. N validation points are 
drawn randomly according to P] ~ pc{p' m ),i = 1, - . • ■ N. The validation sample set is 
formed 


— { (P m i , (f\ (P ml)t02 (P ml))> • • • (Pm.Vi^l ( P m .V ) ■ 02 (Pmjv))}> Pmi ~ P(Pm)- 

(5.18) 

by appealing to the truth outputs at each of the Nc input validation points. The size of 
the sample set is given by the validation sampling theorem 


N c 


In £2 

ln(l - £ 1 ) 


(5.19) 


where \z ] is the smallest integer that is greater than 2 . 
The model prediction error is then computed as 


U c = max 

X va 


f max 


l^(P m! )-0 1 (P m ,)| l4(Pmi) ~ fetPm .i 

5oi(Pmi) 5c>2(Pmi) 


(5.20) 
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Where £<?,( p' m ) : -4 J?+ and ff P2 (p' m ) : -»• JR+ are strictly positive scaling func- 

tions. The notion of the the uncharacterized region is defined as 


l^f(Pm) - 0l(Pm)l \<t>2 (Pm) “ ^(Pra 




max 


501 (Pm) 


502 ( Pm) 


< U c 


(5.21) 


The model prediction error is the maximum, absolute difference between the truth and 
surrogate outputs over the set of randomly drawn, validation input points P mi , i = 1, . . . , N. 
The uncharacterized region, T^, is that subset of the PO model design space for which 
the scaled surrogate error is greater than the model prediction error Uc- The validation 
provides no information as to the location of the uncharacterized region nor does it provide 
any insight as to the magnitude of the errors in that region. 

For the eddy-promoter problem, the scalarization parameter w £ [0, 1] used to find 
the PO manifolds has been used (after being appropriately cleaned up to remove duplicate 
points and non-PO segments) to define the probability density function pc{ Pm)- The 
probability density for the randomly selected values of W ~ f\v{w) is uniform. f\v{w) = 1.0. 
The radius ry in Equation 5.16 has been set to r\- — 0.01. The validation sample size Nc has 
been set to Nc = 21 and. according to the sampling theorem in Equation 5.19. S\ = 0.0639 
and S '2 = 0.2500. The scaling functions have been set to (je = 1.0 and g v = 1.0. The 
validation sample set has been formed through appeals to the simulation at each validation 
point and the model prediction error Uc has been computed and found to be Uc = 0.02506. 

It is immediately obvious that the model prediction error found in this section is sig- 
nificantly smaller than the global validation error computed in Section 5.2.1. The primary 
reason for the improvement is that the PO input manifold corresponds (through good for- 
tune) to regions in the model design space fl m for which the surrogates do a very good job 
of approximating the truth. There is no reason to expect an improvement in general, and 
in fact, the opposite would more likely be true. The optimization process has the tendency 
to find the regions of the design space for which the surrogate approximations are the most 
poor. An example of this would be regions in which the surrogates have large undershoots 
that are not physically justified. The optimization process would seek these regions out 
and the surrogates would falsely indicate that they are PO. Therefore, the reduced error 
for the PO validation of this section is mostly attributable to good fortune and should not 
be expected in practice. The burden to construct accurate surrogates over the entire model 
design space f2 m remains. 


5.3 Surrogate Based Design 

The goal of the design problem is to achieve lower values for the two performance metrics 
$1 = < M<Mp' m ),<Mp' m ):P') : A' -* R and $ 2 = $2(^1 (p' m ), <Mp' m )> p') : fl 1 -+ 1R. In 
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most design problems of relevance however, the metrics represent competing goals and can 
not be simultaneously reduced to a single, optimal point. The purpose for pursuing the 
Pareto optimization approach has been to maintain the maximum amount of flexibility as 
to how to competing metrics will be prioritized during the design stage. The Pareto-optimal 
solutions were shown to satisfy a wide class of multicriteria optimization formulations in 
Section 2.4.2. Even with the design space reduced from f l to C A , the solution to most 
problems of interest are still available and any of the formulations given in Section 2.4.2 
apply. Furthermore, the design problem, for the two performance metric case examined 
here, is trivial to analyze as the Pareto analysis reduces it to a single tradeoff curve. 

For the general, two performance metric problem, a constrained optimization problem 
formulation is chosen. The optimal input points p 7 , q G Q = (1, . . . - Q ). can be expressed 
as 


p 7 = arg __ min _ ($i(0i(Pm)- 02(Pm)- P*))- V 7 e Q. (5.22) 

{p'€ni4>.>(oi(p^),Oj(p' m ).p')<'J>''»} 

where is the upper bound for the constraint metric. ^ 2(^1 (p' m )- 02 (Pm)- P*)- T * ie 0 Pt'' 
mizers. p 7 . are easily obtained from the PO output manifold and the impact of a changing 
constraint bound. is easily assessed. 

For the eddy-promoter problem, the formulation that will be pursued will be to find 
the Q designs with the best pressure performance and that achieve a prescribed tempera- 
ture performance, 9 q . q G Q — (1 Q)- The constrained optimization problem can be 

expressed as 

p 7 = arg min _ 'f(^o(Pm)- P)- V</ G Q. (5.23) 

{pen | 0(0~ o (pm).p)<0''} 

where the different temperature performance bounds 6 q , q 6 Q represent what, in practice 
might, be evolving design goals or a tradeoff analysis. 

Three values {Q = 3) of the temperature performance bgund iii Equation 5.23 have been 
chosen. For the first, the design goal is the find minimum temperature performance design 
and, by inspection of the PO output manifold in Figure 5-1, 0 1 = —2.2626. For the second 
design, the temperature performance constraint has been relaxed and is 9 = -2.0239. 

Finally, for the third design. 9 3 = -1.7106. The optimal design point has been found 
for each by finding the value of the scalarization parameter w that corresponds to each 
design and then appealing to the PO input manifolds. The surrogate-predicted, optimal 
input points are p ^ — (0. 1008. 0.2875, 1 .0000, 1 .0000). p - = (0.8080, 0. 7 286. 0.5448. 0.1868). 
and p 3 = (0.8549,0.4729,0.6449,0.0000) and the corresponding surrogate, performance 
metric pairs are flfp 1 ) = (—2.2626,12.5137), fl(p 2 ) = (—2.0239,10.6495), and ^(p , ) = 
(-1.7106,9.4801). The surrogate, performance metric pair at an input point p' is defined 
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(b) Design 2; p 2 = (0,8080.0.7286.0.5448.0.1868). 
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(c) Design 3: p 3 = (0.8549.0.4729.0.6449.0.0000). 

Figure 5-4: Eddy-promoter configurations of the surrogate-predicted. Pareto-optimal de- 
signs used in the design study. 


as 

n(p') = (0(p')-*(p')) = (e(0o(p:j.p^*(^o(p;j,p')) • (5.24) 

The configuration for each of the surrogate-predicted optimal designs is plotted in Figure 
5-4. From the results shown, it is evident that the topology-change capability does impact 
the results. Depending on the preferences given to each performance metric, the optimal 
geometry may consist of a single eddy-promoter inclusion as in Design 1, or two distinct 
bodies as is the case for Designs 2 and 3. However, the results presented here are not 
sufficient to assess whether the two-topology capability enables a significant improvement 
in the optimal results. One way to determine the extent of the two-topology impact is to 
pursue the Pareto analysis over the subset of the design space that corresponds to only 
single-body geometries. By overlaying the single-body, restricted PO output manifold on 
the PO output manifold plotted in Figure 5-1. the impact that the two-topology capability 
has could be quickly assessed. 
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5.4 A Posteriori Error Analysis 

The a posteriori error analysis step of the surrogate-Pareto framework is identical to that 
of the baseline framework described in Section 4.4. but for the predictability analyses, 
the design space is restricted to narrow “tube" P. c encompassing the surrogate PO input 
manifold C A . Two error estimate formulations for the predictability analysis are examined 
in this section: the proximal region analysis in which error bounds are valid for a finite 
sized region near the surrogate-predicted optimizers and the proximal candidate analysis 
in which the bounds are valid for a specific, randomly chosen design near the surrogate- 
predicted optimizer. The optimality analysis proceeds very similarly to that of the the 
baseline framework describe in Section 4.4.3 but is slightly modified to take into account 
the full Pareto family of optimal designs. 

The three analyses, proximal region, proximal candidate, and optimality are examined 
in the following three sections. The format for each section is to present the theory for the 
general, two output, two performance metric problem and then to give the corresponding 
results for the two-bodv. eddy-promoter example. 

5.4.1 Predictability — Proximal Region 

To begin the predictability analysis, it is first assumed that a prediction region V' x C P L 
has been defined that contains the surrogate-predicted optimizer p*. As was done for 
the baseline analysis, the surrogate-predicted optimal input design vector p* has the form 
p* = (p* n .p*) where p>‘ m is the vector of modeled inputs and p* is the vector of inputs for 
which the performance metric response is known analytically. The prediction region can 
similarly be expressed as Vt x = x p*. The p-measure of the prediction region relative 

to the PO design space is computed as 

HpA V m Sl ) = / _ Pc(Pm)dp m = Si, (5-25) 

JV m' l 

and. for the analysis to be valid, must be at least as large as An example of how to 
precisely define and measure the prediction neighborhood is given for the eddy-promoter 
problem later in this section. 

With the prediction neighborhood defined, probabilistic bounds on the truth output 
values for inputs points that are in a finite-sized, subset of V* x are developed. The resulting 
predictability statement for the truth outputs follows: With probability of at least 1 - £ 2 , 
there will exist points p 1 6 f C Vt x such that 


Id: < <Pl(p'm) < “op 

( 5 . 26 ) 

I 02 < AiiP'm) < u oi- 

( 5 . 27 ) 
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where the truth-output bounds are based on the validation results and are 


= min [<MPm) -Ucg^iPm)], (5-28) 

{P } 

= max .[0l(Pm) + 0/:$0i(Pm)]* (5-29) 

{p'SPm-j} 

4a = , min [02 (Pm) -^Sda(Pm)]. (5.30) 

{ P'6Pm fl } 

u 02 = , max f^2(Pm) + ^£5d 2 (Pm)] ! (5.31) 


and (7^ is the PO model prediction error computed from Equation 5.20. The bounds given 
in Equations 5.28 — 5.31 can be computed entirely from the surrogate values, require no 
additional truth calculations, and, therefore, are very inexpensive to evaluate. 

Equations 5.26 — 5.31 effectively bound (within a confidence of 1 — £ 2 ) the surrogate 
output error for points “near ’ the surrogate predicted optimizer in terms of the outputs. For 
the bounds to be useful to the designer however, they need to apply to the truth performance 
metrics. $l {<£[ (p 7 m ). ^(Pm)? P ; ) anc ^ ^(^[(Pm)* - PO The predictability statement 

on the performance metrics reads; With probability of at least 1 — £ 2 - there will exist a 
region of non-zero measure F' C in the neighborhood of p* such that for all points 

p' 6 n 

£*l < *l(^r(Pm).^2(Pm)-P') < U*,. (5.32) 

< <M<tf(pU^(p'J-P') < U* 2 . (5.33) 

where the performance metric bounds are 

£$, = „ min , $i(zi,Z2-p'). (5.34) 

{p^v^^eZuzi&Zi) 

U<t>. — max $ 1 ( 21 , Z 2 . p'), (5.35) 

(P'€'P; 1 ,2i€2i 1 2262 2 } 

L< j>., = min ^* 2 (^ 1 , 22 , p')- (5.36) 

- {p'€^;,,*i €2 m*€2j} 

= max $ 2 (zi, Zo. p'). (5.37) 

{p'e'P: 1 'Zi€Z l ,z i eZ2} 

and the output ranges, Z\ and are determined by Equations 5.28 — -5.31 and are 

Z x = {z\l 0l <z<u 9l }, (5.38) 

Z 2 = {z \ Ifr < z < Up 2 }. (5.39) 

As was the case for the output bounds, the truth performance metric bounds given in 
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Equations 5.34—5.37 require no additional truth calculations, and therefore, are inexpensive 
to evaluate. 

To define the prediction neighborhood V‘ x precisely, the distance metric A(p ml .p m2 ) 
used in Section 4.4.1 and given in Equation 4.31 has been used. With this definition for 
the metric, the prediction neighborhood around the surrogate predicted optimal point p q m 
is then defined as that region 77, of size £\. that minimizes 


r A = max A(p m .p£J. 

P 


(5.40) 


For the eddy promoter problem, the distance metric function has been defined based on the 
outputs and is selected from the general form 


A(p ml .p m2 ) = max 


l^o(Pml) - 0o(Pm 2 )l l^o(Pml) ~ <A)(Pm 2 ) 


hg 


K 


(5.41) 


where kg and h v are positive scalars. The pressure output scalar h v in Equation 5.41 has 
been set sufficiently high so that the the neighborhood with minimum 9 0 sensitivity results. 
The distance metric in this case reduces to 


A(p ml - Pm 2 ) = |0o(Pml) -0o(Pm 2 )|- (5-42) 

The output bounds in Equations 5.2S — 5.31 have been evaluated for each of the three 
designs found in Section 5.3. The output bounds have then been used to evaluate the bound 
quantities given in Equations 5.34 — 5.37. The resulting predictability statements for the 
truth values of the performance metrics are given below. For Design 1. the statement reads: 
With probability of at least 75%, there will exist points in T 1 C Vl x such that for p' € T 1 


12.1440 < Vtyl (p' m ). p') < 12.5841 (5.43) 

-2.2816 < Q(dJ (p' m ), p') < -2.2359. (5.44) 

For Design 2. the statement reads: With probability of at least 75%, there will exist points 
in F 2 C V'l x such that for p' G T 2 

10.5211 < V(ipo(p'rn)’P') < 10.7510, (5.45) 

-2.0631 < 0(^(p' m ).p') < -1.9880. (5.46) 

Finally, for Design 3. the statement reads: With probability of at least 75%. there will exist 
points in T 3 C such that for p' 6 T 3 

9.1360 < (p' m ), p') < 9.6281. (5.47) 
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Figure 5-5: Output plane plot of the three design point, surrogate-performance pairs. 
ri(p'), i = 1.2.3 (•), and the proximal region predictability boxes for each design. 


-1.8010 < 0(0 o r (p!„)-p') < -1-6361. (5.48) 

The predictability statements are valid for each design for the corresponding prediction 
neighborhood. 

The results given above have been plotted in the output plane in Figure 5-5. In the 
figure, the surrogate performance metric values for each design are plotted as solid dots, 
the bounds on the truth performance metrics are plotted as boxes, and the surrogate PO 
output manifolds for the eddy-promoter and the plane-channel heat exchanger are shown as 
solid and dashed lines , respectively. This figure can be compared to the baseline surrogate 
framework results given in Figure 4-3 in which the boxes were significantly larger. The 
bounds in Figure 5-5 are quite sharp with a clear discrimination between each of three 
designs in terms of both the pressure and temperature performance metrics. Additionally, 
the surrogate predictability boxes indicate the the surrogates are accurate enough to clearly 
show (with 75% confidence) that the eddy-promoter heat exchanger performance is better 
than the plane channel performance for all three designs. 

Several factors contribute to the improvement in the predictability boxes plotted in 
Figure 5-5 relative to the baseline framework results plotted in Figure 4-3. First, the model 
prediction error was significantly lower for the Pareto validation over f than for the 
global validation over Q m . This is mostly good fortune, as PO input manifold happens to 
lie in regions of for which the output surrogates are accurate. The second reason is 
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due to the increased predictability that can be contributed to reducing the dimension of 
the validation space from three (the three model inputs of the model design space 9., n ) to 
just a single input (the scalarization parameter w). A third and final contributor to the 
sharp predictability boxes is the input p a = for which the performance metric responses 
are known analytically. This effect was present in the baseline results as well, but was 
overwhelmed by the very poor model validation error computed in Section 4.2. 


5.4.2 Surrogate Predictability — Proximal Candidate 


The predictability can be evaluated in terms of a randomly selected design point near the 
surrogate-predicted optimizer in much the same manner as it is for proximal candidate 
analysis of the baseline surrogate framework. The only difference is that the prediction 
neighborhood for the Pareto framework analysis is a subset of one-dimensional PO design 
space 9 C instead of the full design space 9. The remainder of the analysis is identical to 
that presented in Section 4.4.2. 


Following the presentation of Section 4.4.2. a prediction region V * C 9 C of p- measure 
a and that contains the surrogate-predicted optimal point p* is defined in a manner similar 
to that presented for the Proximal Region analysis of Section 5.4.1. A sample candidate 
design is selected randomly according to P*, ~ p Pc ( p m ). The density function pp c (p m ) 
is defined as 


PV c (Pm) 


~P(Pm 

a 


. Vp m €V* ma . 

V 


(5.49) 


where V* na C 9f n ' is the sub-manifold of V* that corresponds to the model design space 
and the full neighborhood is the tensor product V * = { 'P^ na x p‘). The full candidate 
design vector P* is given as P* = (P* n . p* ( ) where p* is the vector of analytic inputs of the 
surrogate predicted optimizer p* = (p,VPa)- 

A second small parameter, e c , is introduced and is related to the p-measure of V '^ na . a. 


by 


= 


.(l_(l_ CT )^-l) 


cr(Nc + 1) 

where 0 < £ c ,cr < 1, and Nc is the validation sample size given in Equation 5.19. 


(5.50) 


Now. given the two inputs to the analysis s c and a, and the model prediction error Uc - 
the following statement can be made: With probability of at least 1 — £ c , the truth output 
values for the input vector P*„ are bounded by 


^ <0[( Pm)<<, (5.51) 

l C o, < < 4- (5-52) 
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where the bound values are 


^ =MK)-Uc9* l (*mh 

(5.53) 

< = wp;)+^,(a 

(5.54) 


(5.55) 

4 = ^(p;) + t/^ 2 (p;). 

(5.56) 


Equations 5.51 — 5.56 bound the truth value of the outputs for the candidate design point 

P* 

x m ' 

The analysis can be extended to the performance metrics the same way as it has been 
in the previous sections. The following predictability statement results: With probability 
of at least 1 — £ c , the truth performance metric values at input point are bounded by 

L%,<*x(<pJ( P m m ),4(P* m ),P’) <U^ (5.57) 

L% 2 < <5,(0[(P^).^(P; n ).P*) < l% 2 . (5.58) 

where the performance metric bounds are 



= min 

{.-i€3f.sa€3S} 


Z? , 

P*)- 

(5.59) 

u% x 

= max 


%2' 

p*) : 

(5.60) 

L%, 

= min 


Z'2 . 

p*). 

(5.61) 


{«i€2f,sa62S} 


ui 2 

= max 

^2(21, 

Z‘2 • 

p*). 

(5.62) 


and the ranges for the outputs, Z\ and Z$, at the candidate design point are obtained from 
Equations 5.53 — 5.56 and are 


2i c = {*|^ <*<*&}, (5.63) 

= {z I l C 6, < 2 < «£,}- (5.64) 

The proximal candidate analysis has been carried out for each of the three eddy- 
promoter. heat exchanger design points selected in Section 5.3. Given the Nc = 21 vali- 
dation points used for the PO design validation, the values of s c = 0.2500 and a = 0.1795 
are obtained from Equation 5.50. The distance metric used to define the neighborhood V* 
is identical to the one used for the proximal region analysis Section 5.4.1 and is given in 
Equation 5.42. Monte-Carlo sampling has been used to obtain the three candidate designs. 
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(a) Candidate Design 1: P 1 = (0.1004.0.3018.0.9950.1.0000). 
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(b) Candidate Design 2: P 2 = (0.8062.0.7290.0.5434.0.1868). 
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(c) Candidate Design 3: P ! = (0.8083.0.7211,0.5527.0.000). 

Figure 5-6: Eddy-promoter candidate configurations for the three surrogate-predicted. 
Pareto -optimal designs used in the design study. 




A single periodicity cell for each of the three candidate designs is given in Figure 5-6. The 
input values for each of the candidate designs are P 1 = (0.1004.0.3018.0.9950.1.0000). 
P 2 = (0.8062.0.7290.0.5434,0.1868). and P 3 = (0.8083.0.7211,0.5527.0.000). . , 

The candidate design geometries plotted in Figure 5-6 can be compared to the surrogate- 
predicted optimal design geometries plotted in Figure 5-4. The candidate and surrogate- 
predicted geometries for Design 1 are very similar, but do differ slightly. For Design 2. the 
candidate and surrogate-predicted geometries are nearly identical. There is a significant 
difference between the candidate and surrogate-predicted geometries for Design 3: in fact, 
the candidate geometry for Design 3 is almost identical to Design 2. The reason for the 
apparent preference toward geometries similar to the surrogate-predicted Design 2 is a 
consequence of using the scalarization parameter to randomize the design space in the 
Pciv'm) definition. Because the model inputs p m remain nearly constant over a large portion 
of the range of w (see Figure 5-3), the effective pc(p' m ) corresponding to this range is very 
large. The Design 2 input point falls into the single geometry range of w and Design 3 is 
near the. Because pciv'm) is ver y lar S e over t ^ 1 ‘ s ran § e - there is a very high probability that 
randomly selecting points near Designs 2 and 3 will result in the candidate designs that we 
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Figure 5-7: Output plane plot of the surrogate response for the three candidate design points 
(o), the surrogate-predicted optimizers (+), and the proximal candidate predictability boxes 
for each design. 


obtained. Design 1 is also very close to the Design 2 geometry in Figure 5-3 but the less 
likely outcome of selecting a candidate similar to Design 1 occurred. 

The results of the analysis are given in terms of the predictability statements. For 
Design 1, the candidate predictability analysis reads: With probability of at least 75%, the 
truth performance bounds for the randomly drawn candidate design point P 1 are 

-2.2805 < 0(^(P^),P l ) < -2.2433. (5.65) 

12.4214 < (P’jii.) < 12.5549. (5.66) 

The surrogate-predicted performance metric values for P 1 are 0* = 0(#o(Pm)-P 2 ) = 

—2.2616 and 'P* = — 12.4933. The bounds are plotted as predictability 

boxes in Figure 5-7, along with the surrogate performance pair outputs at the optimal 
design point nip 1 ) = (Q 1 ,# 1 ) (plus sign) and for the candidate design ft(Pi) = (0 1 ,$ 1 ) 
(open dot). 

For Design 2, the candidate predictability analysis reads: With probability of at least 
75%, the truth performance bounds for the randomly drawn candidate design point P 2 are 

-2.0630 < S(9j {P 2 m ). P 2 ) < -1.9881, (5.67) 
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(5.68) 


10.5215 < <6{^(P 2 ),P 2 m ) < 10.7484. 

The surrogate-predicted performance metric values for P 2 are 0 2 = 0(0o(P^),P 2 ) = 
-2.0240 and T 2 = 'P^ofP^Pm) — 10.6496. These results are plotted in Figure 5-7 as 
well. 

Finally, for Design 3. the candidate predictability analysis result reads: With probability 
of at least 75%, the truth performance bounds for the randomly drawn candidate design 
point P 3 are 

-1.8307 < 0(0j( P^),P 3 ) < -1.7322. (5.69) 

9.4284 < T(t/^(P 3 ).P^) < 9.6568. (5.70) 

The surrogate-predicted performance metric values for P 3 are 0 3 = ©(#o(Pm)- P 3 ) = 
— 1.7786 and T 3 = 'P(i/»o(P 3 ), P'^) — 9.5574. The results are plotted in Figure 5-7. 


5.4.3 Design Optimality 

The global validation results provide the information necessary to estimate the amount, of 
effort required to improve a design beyond a computable amount determined by the the 
global model prediction error. Based on the global model prediction error. L . the maximum 
absolute output error over a 1 — £\ fraction of the design space is known with a confidence 
level 1 — S‘ 2 * For the error at a point to be worse than the model prediction error, the point 
in question would have to lie in the uncharacterized region T. This notion is used in the 
development of the optimality estimates below. 

To find the effort required to improve upon a design, an expanded achievable set A 3 is 
first introduced 


-4° = {s € R 2 
where 


and 


3p' € Q',z\ € Z^ { ,Z 2 € Z$ 2 s.t. <f>i( 2 :i. r- 2 - p') < «i , 4 * 2(21 • < •'>■ 2 }. 



(5.71) 
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lo 1 = . min } 0i( p m) - Ug Pl ( p)- 

(5.74) 

= max 4>i(Pm) + Ugo x (p'). 

(5.75) 

1 £ = , min , faiPm) ~ u 9o *(P')- 
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(5.76) 

= { max ; 02(p' m ) + Ug*(p'). 

(5.77) 
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The boundary of A 0 is denoted dA° and is the expanded, Pareto-Optimal output manifold. 
The interest is in how much effort it would take to find a point p' that has performance 
metric response not in A° . For a point p 1 to exist such that 

(*i(£(Pm), £(Pm). p')> fc(p'J, p')) 0 -4°, (5.78) 

it would have to lie in the uncharacterized region T. 

The remainder of the analysis is similar to the baseline surrogate optimality analysis that 
was presented in Section 4.4.3, A sequence of iV opt , independent, identically distributed, 
randomly selected inputs, ,P mvop{ , is drawn according to the probability density 

function p(p m ): P mi ~ p( p m )»i = 1, . . . , A ropt . The optimality statement then follows that, 
the probability that at least one of the randomly selected design points will have truth 
performance values that lie in A unch is less than 1 — Sl. The number of random input 
vectors drawn, iV opt , is computed as 

N opt < jv (T - l) , ( 5 .79) 

where the parameter ei E]0, 1[ is set independently. The optimality result states that, with 
a confidence of less than 1 — El, N opl additional truth evaluations at randomly selected 
points will produce a design better than the surrogate selected design. This is an estimate 
of how much work would have to be expended to improve the design beyond a given amount. 

The optimality analysis has been performed for the eddy-promoter problem. The PO 
boundary of the expanded achievable set. A°> has been computed and is plotted as a dashed 
line in Figure 5-8 and the surrogate PO output manifold is plotted as a solid line. The 
parameter ei in Equation 5.79 is set to ei = 0.333 and, given the N = 24 validation points 
that have been sampled in Section 5-8, the value for for iV opt is jV opt = 48. The optimality 
statement follows: If an additional 48 input points are drawn randomly according to p( p), 
P m t ~ p{Pm ) • i = 1,...,48, the probability that at least one will have truth output pair 
n r (P m .) = ('I'(^(P m ,),P 1 ),0(^(P m ,),P 1 ) 0 A° is less than 66.7%. 

What the optimality statement gives is an estimate to how much additional work (48 
additional truth evaluations in this case) would be required to have at least one point that 
achieves a performance better than A° with confidence of 1 — sl. For the eddy-promoter 
example, the distance between dA and <9*4° is very large, suggesting that the potential 
improvement is likely worth the additional computational effort. Secondarily, the wide gap 
between <9.4 and <9-4° indicates that the output surrogates are not performing well enough 
to provide meaningful results in terms of global optimality. Although the predictability 
for the family of points selected by the Pareto optimization were acceptable according the 
predictability analyses in Sections 5.4.1 and 5.4.2, the surrogate predicted PO family of 
input points are not likely close (in performance) to the truth PO family. 
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Figure 5-8: Output plane plot of optimality results for the surrogate Pareto framework. 


For a final estimate of the surrogate performance, a hypothetical optimality analysis can 
be performed to determine the magnitude of the global model predication error U that would 
be necessary to achieve acceptable results. This has been performed for the eddv-promoter 
problem by assuming a model prediction error Uh yp that is lower than U = 0.0968 by a 
factor of 5. U h , )P = U/o = 0.0968 = 0.0193. With the hypothetical model prediction error 
U^p = 0.0193. the associated expanded achievable set A° hyp has been defined as in Equation 
5.71 and the expanded PO output manifold 0A° hyp has been determined. The results from 
this analysis are plotted in Figure 5-9. In Figure 5-9. true expanded PO output manifold. 
0A°. is plotted as a dashed line, the hypothetical expanded PO output manifold. dA° hyp . 
is plotted as a dashed-dotted line, and the surrogate PO output manifold DA is plotted 
as a solid line. The very narrow gap between dA°h yp and dA suggests that if the output 
surrogates were improved by a factor of 5, the surrogate designs may be close enough to 
the truth optimal designs to be acceptable. 

The optimality analysis highlights how poorly the surrogates approximate the truth 
outputs for the eddy-promoter. The construction point density should be sufficient to 
approximate relatively smooth, well-behaved outputs; N co = 256 translates roughly to 6 
points in each of the 3 input directions for a full factorial design. However, the truth out- 
puts are far from smooth for the eddy-promoter: in fact, for changes in the model input 
variables, the outputs undergo subcritical bifurcations in which the output is. possibly, dis- 
continuous. The discontinuities in the input-output function are very difficult to precisely 
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Figure 5-9: Output plane plot of hypothetical optimality results (dotted line) for a factor of 
five improvement in the model prediction error, the actual optimality bound (dashed line) 
and the surrogate Pareto output manifold (solid line). 


isolate in three input dimensions and produce large surrogate dispersions from the truth 
output. No effort has been made in the global surrogate construction step to account for the 
possibility of discontinuous outputs and the radial basis interpolation method exacerbates 
the difficulties in these regions. A significant improvement in the surrogates (comparable 
to what was assumed for the hypothetical analysis above) would likely require that discon- 
tinuity locations be determined precisely, and accounted for separately. One approach is 
use linear interpolation across the discontinuities as has been done for the lift coefficient 
surrogate described in Section 6.4.1. 
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Chapter 6 


Experimental Airfoil Optimization 


In this chapter, the application of the baseline surrogate framework to the experimental 
optimization of a multielement airfoil is presented [54]. Although the application discussed 
in this chapter is experimental and. therefore, has a zero-mean, measurement error associ- 
ated with the outputs, the errors are assumed to be small and the outputs are treated as 
noise-free. The surrogate framework has been extended to correctly handle noisy outputs 
[78]. but these techniques are not used here. 

The model has been developed, and the data generously provided, by Landman and 
Britcher at Old Dominion University (ODU) [44. 45]. The motivation for the development 
of the experimental model central to this chapter is the inherent difficulty in examining many 
design points experimentally. The model is a three-element airfoil model with internally 
embedded actuators (Figure 6-1). It has a nested chord of c = 18 in., a span of b = 36 
in., and was designed for low-speed testing in several wind tunnels, including the NASA 
Langley Research Center 2- by 4-foot tunnel and the ODU 3- by 4-foot low-speed facilities. 
The main element chord is c mam = 14.95 in., and the flap and slat chords (expressed as a 
percentage of the nested chord) are 30 and 14.5 percent, respectively. The flap and slat are 
both deflected to 30° for all tests. The particular model used for the test reported here has 
been developed for low Reynolds number testing to prove the experimental testing concepts. 
However, the techniques that have been developed should be applicable to higher Reynolds 
number testing as well. 

The flap actuators are computer controlled and position the flap horizontally and ver- 
tically (x and y, respectively). The model has been used in the ODU tunnel to compile 
baseline values for lift coefficient C; versus flap gap and overhang at fixed angles of attack 
and slat riggings. A first-order optimizer that uses a variant of the method of steepest 
descent [24. 10] has been demonstrated in real time[45]. 

The capability of the computer controller to automatically take data at a prescribed set 
of ( x.y ) trailing edge flap coordinates makes the airfoil model setup ideal for application 
of the surrogate methods to an experimental problem. Additionally, the opportunity to 
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Figure 6-1: Three-element model with internal flap actuators. 

compare the surrogate results to those obtained with in on-line optimizer is unique to this 
experimental test. 

The baseline, nonparametric-validated surrogate framework presented in Chapter 4 is 
applied to the airfoil optimization problem. The problem examined here is a single output 
problem that is a function of only two inputs; the (,r. y) coordinates of the trailing edge flap. 
The baseline surrogate framework provides a practical means to incorporate experimental 
data directly into the design optimization process. The off-line surrogate approach to 
the design optimization has several advantages to on-line optimization [48. 79]. First, by 
construction, surrogates are computationally inexpensive and are thus easily incorporated 
into optimization procedures. Additionally, the low computational requirements create a 
highly interactive and flexible design environment, which allows the designer to easily pursue 
and examine multiple design points. Second, the number of appeals to the experiment or 
simulation is known a priori , which ensures that the design can be accomplished without 
exhausting available resources. Third the surrogate approach offers a natural means to 
incorporate data from previous runs and/or other sources. 

As regards disadvantages, the primary drawback is that in high dimensional design 
spaces, surrogate construction is difficult and design localization is poor. The poor design 
localization has been demonstrated for the baseline surrogate framework optimization of the 
eddy-promoter heat exchanger in Chapter 4. A second limiting factor in the application of 
the surrogate approach to experimental tests is the need to validate the surrogate at input 
points chosen randomly in the design space. This capability, present in the experiment 
central to this work, is not typical of most experimental tests. Finally, surrogate-based 
optimization introduces a new source of error. The surrogate validation strategy and error 
norms discussed in this paper seek to quantify the discrepancy between the surrogate and 
the experiment by providing estimates to the system predictability and optimality. 
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Figure 6-2: Definition of gap and overhang. 


In this Chapter, the experimental model and testing methods used are briefly described 
first. Second, the optimization problem that is central to the work is presented. Third, 
the three steps steps of the baseline surrogate framework (i.e.. construction/ validation, 
surrogate-based optimization, and a posteriori error analysis) are presented, the inputs 
to the framework are summarized, and an overview of the more sophisticated surrogate 
algorithms is provided. Finally, sample results obtained from the surrogate framework for 
output, maximization and multiple-target designs are presented, and a comparison between 
the surrogate approach and direct insertion results reported previouslv[45j is provided. 

6.1 Experimental Testing Methods 

An important practical problem encountered in wind-tunnel testing of multielement airfoils 
is the need to test a range of configurations to ensure that the optimum is selected. Un- 
fortunately. this testing can be prohibitively time consuming if one considers all possible 
variables, such as flap position and deflection, slat position and deflection, overall angle 
of attack, and Reynolds number. For example, a range of flap locations and orientations 
relative to the main element is typically tested. In a cryogenic or pressurized facility, model 
geometry changes necessitate large delays in testing. These delays often result in investi- 
gators choosing a sparse test matrix and an optimum that is based on only a few points. 
The ability to move the flap under computer control provides a unique opportunity to ex- 
plore the entire range of useful gap and overhang values (Figure 6-2). Two typical pressure 
distributions are shown in Figure 6-3. where the ordinate is the pressure coefficient C p and 
the abscissa is distance from the leading edge expressed as a percent of the nested chord. 
The data for Figure 6-3(a) represents a point near the peak C; for this configuration, and 
the plot in Figure 6-3(b) indicates full separation over the flap. 

In this experiment (performed by Landman and Britcher). the flap actuators, tunnel flow 
setting, and data acquisition were controlled by a personal computer running Lab View [50] 
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(a) Fully attached flow over all elements. 


(b) Detached flow over the flap. 


Figure 6-3: Experimental pressure data. 


software. A program was written to allow any number of flap positions (in x and y) to 
be sampled in any order. Wind tunnel power was controlled such that at the beginning 
of each test the tunnel was restarted to avoid hysteresis effects [45]. The experimental 
setup allowed the user to start the program, which at each location in turn automatically 
measured the free-stream properties, sampled and recorded pressures around the centerline 
of the model, and then calculated lift coefficients for the three-element airfoil. This process 
required approximately 2 min. for each data point. 

Test matrices were developed to survey flap positions, which ranged from approximately 
0.8 - 3.5 percent (gap) and —0.4 - 3.4 percent (overhang) relative to the nested chord c. Two 
angles of attack and two slat geometries were selected. An angle of attack a of 8° was chosen 
as representative of an approach value. An a of 14° represented the limit of good-quality 
two-dimensional flow for the ODU tunnel installation without tunnel wall boundary-layer 
control. Two slat settings were chosen: a slat gap of 3.03 percent with an overhang of 2.46 
percent and. for a smaller gap setting, a slat gap of 2.17 percent with a slat overhang of 
— 1.46 percent. 

Positional accuracy was enhanced by requiring that the flap move to a reference point 
above and behind the desired evaluation points ( x re j > x eva i . y re j > y eva i ) and then 
back to the evaluation point. This eliminated any effect of backlash in the mechanical 
drive-train. Two simple tests provided an indication of the inherent collective error due 
to instrumentation and positioning. The first test involved two separate evaluation points: 
the first point was in a region in which the flow was known to be fully attached to all 
elements, and the second point was chosen in a region in which flow over the flap was fully 
separated. The positioning program was used to move the flap between a reference point 
and one of the evaluation points. The tunnel was restarted before every evaluation, and the 
test was repeated 30 times in each case. The standard deviation of C[ was found to be 0.004 
for the separated case (0.16 percent) and 0.0118 for the attached case (.36 percent). For 
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the second test, the program automatically sampled 29 points over the entire test region 
for two different trials. The error in C{ between the two runs averaged 0.71 percent with 
a standard deviation of 0.75 percent. Although these tests are not exhaustive, they do 
provide a benchmark for the C/ error. 

The turbulence intensity in the ODU tunnel was measured at less than 0.2 percent. 
Flow quality over the model was monitored through 12 spanwise taps: 6 on the flap, and 
6 on the main element. The flow was considered to be two-dimensional if the magnitude 
of the spanwise nonuniformity was less than 5 percent of the total C p variation over the 
entire model [49]. The data presented are uncorrected for boundary effects were taken at a 
Reynolds Re number of 1 x 10 6 based on the nested chord. 

6*2 Optimization Problem 

We begin by introducing a vector p of M design inputs that lie in the input (or "design*) 
domain Q C JR‘ U , an input-output function S(p) : 0 -> iR, and an objective function 
T(S(p),p.A) that characterizes our design goals, where A is a vector (or possibly scalar) 
design parameter. For the work presented here, we set p = {x.y) (the x- and //-positions 
of the flap) as the M = 2 inputs and restrict ourselves to an input domain Q of reasonable 
flap positions (described in more detail in the results section). The output of interest is the 
lift coefficient. «S(p) = Ci(x.jj). The objective function is T(5(p).p. A) = |S(p) - A| which 
has been referred to as the "discrimination" problem [62]. 

With the above terms defined, the minimizer p* = (x*,y*) to the exact optimization 
problem is given by 

p* = arg min |«S(p) - A| . (6.1) 

In this formulation, the goal is to find that (or "an*') input vector p* = ( x *. y *) that achieves 
as closely as possible the target lift coefficient value A. If the target lift coefficient A is set 
sufficiently small (large), the formulation describes the output minimization (maximization) 
problem, assuming that S{ p) is bounded from below (above). 

In the on-line approach, the experiment is invoked at every optimization step needed to 
solve Equation (6.1). In the off-line approach, a surrogate. <S(p) ^ S(p), for the experiment 
is inserted into the optimization problem. The minimizer. p* = (x*,y*), for the resulting, 
surrogate-based, discrimination problem is then given by 

p* = argmin|«S(p) — A| . (6.2) 

pen 

Here, the optimization proceeds exactly as it would for the on-line approach, but the lift 
coefficient surrogate 5(p) is invoked instead of the experiment. The surrogate problem that 
corresponds to Equation (6.1), but with a general objective function T(S(p). p. A), has been 
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reported by Ye^ilyurt [77] and Yegilyurt and Patera [79]. 


6.3 Surrogate Framework 

The surrogate approach is broken into four steps. In the first stage, surrogate construction, 
experimental results and/or prior information are used to construct the approximation, 
S{ p) «5(p). In the second step, surrogate validation, additional queries to the experiment 
are used to validate the approximation. In the third step of the process, surrogate-based 
optimization, solutions to surrogate optimization problem of Equation (6.2) are obtained. 
In the fourth and final step, a posteriori error analysis, the results of the validation are 
used to analyze the consequences of the surrogate-for-simulation substitution. In the fol- 
lowing subsections, the four steps of the baseline surrogate framework are presented and 
the designer inputs to the framework are summarized. 

6.3.1 Surrogate Construction 

A lift coefficient surrogate S( p) = A{X co ) % 5(p) is constructed using an approximation 
scheme. A : 1R)^ C ° — > L oc (f2) and a construction sample set of input-output pairs 

= {( Pi ,Jl Pl M = 1 N co ], (6.3) 

where R Pi — (?/(£*, y*) is a realization of the experimentally measured lift coefficient for 
the input flap position p ? = and N co is the number of input-output pairs in the 

construction sample. Although the general surrogate framework can handle noisy outputs 
[78]. the noise contribution is neglected in the work presented here. Information from prior 
studies, outside sources, or asymptotic behavior can also be incorporated into the approx- 
imation process either through the definition of the approximation scheme or by including 
the input-output pairs in the construction sample. As already stated, the surrogate frame- 
work makes no assumptions in regard to the approximation technique and will accept, 
and assess, any approximation A{X co ). Also, no restriction is placed on either N c0 or the 
distribution of the construction sample. 

6.3.2 Surrogate Validation 

To proceed with the description of the surrogate validation, and importance function p(p) 
is first introduced. The importance function serves as a probability density function for the 
selection of the validation points: 

/ P(p)dp = 1- (6.4) 

j n 
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The importance function also leads to the notion of a p-measure associated with p(p): for 
any subdomain V C fi, 

^i p {V) = J^p{p)dp K (6.5) 

The p-measure of V is simply the weighted relative M-volume of V. 


With the importance function p(p) defined, the validation sample set is formed 

X w = {(Pi,i?p t ),i = P,~p(p), (6.6) 


where the input flap positions P, for the validation sample set arc drawn randomly according 
to the probability density function p(p). In Equation (6.6), the ~ is as ‘‘is drawn according 
to the probability density function.” The validation sample size N va is given by 


N va 


In c '2 

In ( 1 - £i)' 


(6.7) 


and Ei and e 2 are the two uncertainty parameters described below. The model prediction 
error U is computed from the validation sample set X va as 


U = max 

P.G.V- 


\Rp,-S{Pj) 
9( Pi) 


( 6 . 8 ) 


where Tp" denotes the input points of the validation sample set and p(p) is a strictly 
positive, error-scaling function described in more detail below. 


The result of the construction/validation process is a probabilistic statement that de- 
scribes the global quality of the surrogate «S(p). The validation statement can be compactly 
written as 

Pr{fi p (T) <Ei}>l-e 2 . (6.9) 

where Pr{event} is the "probability of event" and T C fi is the uncharacterized region 
defined as 

T = {p € n\ |<S(p) - <S(p)| > Ug( p)} . (6.10) 

The p-measure of the uncharacterized region is bounded by e 1: and the significance level of 
the nonparametric statistical bound is e 2 - This result can be readily proved [77] with order 
statistics [15] and is included in Appendix B.l. 


For the simple case of g{p) — 1, Equation (6.9) states that, with probability greater 
than or equal to 1 - e 2 - the surrogate error is bounded by U over a region of fi of p-measure 
greater than 1 - E\. Although this statement is suggestive, it gives neither an indication as 
to the location of T nor the magnitude of the surrogate error in T. 
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6.3.3 Surrogate-Based Optimization 

For the optimization problem, it is assume that the design is given Q target drag coefficient 
values A 9 , q E Q = The goal of the optimization is to find the surrogate- 

predicted flap positions that minimize the objective function. 

P q = p*(A 9 ) = argmin|S(p) - A’l, 'iq E Q. (6.11) 

pen 

The Q targets could represent different target lift coefficients during the flap deployment 
schedule, or reflect the goals at different flight conditions (e.g.. take-off and landing). 

6.3.4 A posteriori Error Analysis 

To present the predictability results, the notion of a prediction neighborhood is first in- 
troduced. A distance metric A(a. b) is defined for all (a.b) € ft x ft. which determines a 
"distance" between two input points a and b. Then for any subdomain Pcfi the radius 
of V about a point p is defined as rp(p) = max p / G x> A(p. p'). The prediction neighborhood 
located at point p with a p-measure of 2 . P(p, z ). is that (or a) region V C ft of p-measure z 
that minimizes rp( p). It is assumed that p lies inside 'P(p.z) and that 'P(p.zi) C P(p, Z 2 ) 
for Z\ < Z). It can then be stated that, with probability greater than 1 — £ 2 . for all q £ Q. 
regions F* 7 C V(p q .z[) of nonzero measure exist such that for all p' £ T q . 

|5(p')-5(p’)|< e (p^. (6.12) 

It now remains to bound e(p (? ) and make precise the extent of T q . 

Several bounds are possible on e(p 9 ). which is denoted the predictability gap. The 
predictability of each design can be bounded individually, to obtain 

e(p 9 ) < £(p q : £i), V<7 € Q, (6.13) 

where, for p £ ft and 0 < z < 1. 

£(p.z) = Ug{p,z) + 6{p.z). (6.14) 

and 

g(p.z) = max g(p'), (6.15) 

P'e-P(p.r) 

<J(p,z)= max \S{p') -<S(p)|, (6.16) 

P’SP(p,z) 

and U is the model prediction error from the validation step. Equation (6.8). 

In addition to the joint estimates to the bound on e(p <? ), the the average error over 
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the Q target designs can be bounded. In particular, if it is assumed that the V{p q : £\) are 
mutually disjoint, it can be shown that 


5 £e(p«)<px 


^ n— 1 


where /? = {@\ , . . . , /3q}. and 

L 

C L = {£' G iR^O < /?; < 1, l = 1, . . . , L; £ = 1}, 

/=1 


(6.17) 


(6.18) 


is the set of convex L-tuples. The u nonparametric average" is relevant to multiple-target 
designs and represents the average, as opposed to the worst-case, estimate of the predictabil- 
ity. Also, it is important to note that this predictability bound is calculated entirely in terms 
of the inexpensive surrogate. S( p). 

Finally, for a successful validation (i.e.. n p {T) < e i). the expectation of the size of r 9 can 
be bounded with respect to the validation sample joint probability density. The resulting 
bound is. Vq G Q, 

£ (^| M T ,<, 1 )< 1 + n i _ +[ r i 2 7 - y , ( 6 . 19 ) 

The expression in Equation (6.19) bounds the average /^-measure of the region T 9 , with 
respect to sq. for many validations. 

Several advantages to bounding the errors only to within a finite uncertainty exist [57]. 
First, a sense of stability is obtained in that the estimates apply not only to a single point, 
but to regions F 9 of nonzero measure, assuring that many input points p 9 exist that satisfy 
the error estimates. Second, for the multiple-target case the estimates become sharper 
because there is only a single uncharacterized volume of measure ci- Equation (6.17) is the 
upper bound for the distribution of the single £q-sized uncharacterized region among the Q 
designs. This analysis results in a bound on the average error which is less than the average 
of the individual predictability gap bounds <f(p 9 ,ei). Finally, because the predictability 
analysis is not premised on any particular set of points, the designer has flexibility in the 
choice of the metric A(a. b) (discussed further in the next section). 


6.3.5 Summary of Surrogate Inputs 

To summarize the surrogate framework description, and to highlight the flexibility of the 
environment, we note that four inputs to the process are determined by the user. These are 
listed below: 

i. An importance function p(p) : t 1R+. 
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ii. An error-scaling function £(p) : Cl — > IR+ . 

iii. Two uncertainty parameters, e i and £ 2 , that satisfy 0 < £i ,£2 < 1- 

iv. A pseudometric A(a, b). 

Each input provides the designer with flexibility, and allows the designer’s experience to 
impact and improve the final surrogate-predicted designs. Although poor choices for the 
inputs do not influence the validity of the surrogate results, they greatly reduce the sharpness 
of the results. A short description and explanation of each input follows. 

The importance function p( p) reflects the designers prejudices in regard to the regions 
of that are more likely to contain optimizers. In this context. p( p) is essentially a “prior" 
on p\ To serve this purpose, p(p) is used as the probability density function in the random 
selection of validation points in Equation (6.6). A judicious choice of p(p) (one that is large 
in the regions of the final designs and small elsewhere) can significantly increase the sharp- 
ness of the a posteriori error bounds. The increased sharpness is a consequence of much 
better physical localization (in terms of input variable extent) of the prediction neighbor- 
hood P(p\ci). which in turn reduces the surrogate sensitivity contribution 5(p*.S\) to the 
error bound in Equation (6.14). 

The error-scaling function p(p) can be used by the designer to reduce the impact of 
localized surrogate errors on the error bounds of the final design. Because the model 
prediction error U in Equation (6.8) is global, a large value of p(p) in regions for which 
the approximation is poor will result in a reduced value of the first term on the right-hand 
side of Equation (6.14). provided that the final design does not lie in a region where p(p) 
is large. 

The uncertainty parameters e \ and £2 are related to the number of validation points 
through Equation (6.7). This formula allows the precise budgeting of resources and ensures 
that useful solutions can be obtained. In effect. Equations (6.7)— (6.10) describe what is 
known in a continuous sense about a function based on discrete sampling. Analysis of 
Equation (6.7) shows that, asymptotically for small S\ and £ 2 - N va increases linearly as e\ 
decreases and only logarithmically as £2 decreases. This relationship suggests that although 
we can easily (in terms of validation sample size) increase our confidence in the results 
(smaller £ 2 )? refining the localization of our results (through smaller £ 1 ) is much more 
difficult. The localization has a direct impact on the final error analysis through 5(p,£\) 
in Equation (6.14). The relative difficulty in further refining the localization illustrates the 
need to intelligently select p(p) and where appropriate, A(a, b), both of which can have 
similar effects on the localization error. 

The final input to the surrogate approach is the pseudometric A(a, b). Because A(a. b) 
can be chosen post-validation, various metrics can be examined, and the most appropriate 
selected. One possible trade-off is between design localization (in terms of input variable 
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extent) and predictability in terms of tf(p. sq ) in Equation (6.14). An example of the extreme 
of this trade-off is the sensitivity minimizing metric 


A(a.b) = |<S(a) — 5(b)| (6.20) 

used for the single-point design study of the results section. This metric gives the lowest 
possible d(p.ci). 

6.4 Results 

To demonstrate the surrogate framework, it has been applied to the experimental design 
of multielement airfoils: specifically, the interest is in in the determination of the optimal 
location for the trailing edge flap, based on the lift coefficient C/ in low-speed, high-lift flight 
regimes. The M = 2 design inputs to the problem p = ( x.y ) are the x and y positions 
of the flap, measured from the leading edge of the main airfoil element and normalized by 
the main element chord c mam = 14.95 in. The output of interest is C/. In addition, several 
other configuration and flow condition parameters are fixed for the study. These parameters 
are listed in Table 6.1 and are the Reynolds number Re. the airfoil angle of attack a. the 
flap and slat deflection angles 6f lap and S slal . respectively, and the gap and overhang of the 
slat (expressed as a percentage of the nested chord c — 18.0 in.). 

In this section, the method used for the surrogate construction and the validation re- 
sults are reported first. Second, the single-point design problem of output maximization is 
presented. Third, multiple-target design study is pursued that demonstrates the increased 
sharpness of the nonparametric average error results. Finally, the results of on-line opti- 
mization studies are given and are compared to the off-line, surrogate results. 

6.4.1 Surrogate Construction/Validation 

The construction sample set X c0 consists of 119 input-output pairs that are uniformly 
spaced on a 17 x 7 grid. The {x.y) flap positions for the construction sample are plotted as 


Re 

1,000.000 

a 

14° 

^ flap 

30° 

$slat 

-30° 

9 a Pslat 

2.17% 

overhang s [ at 

-1.46% 


Table 6.1: Fixed design study parameters. 
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Figure 6-4: Surrogate construction points and the input (“design") domain. 

circles in Figure 6-4. The input domain is divided into three subdomains. O = U 
based on the flow conditions over the flap. In the first subdomain f>i, the flow over the 
flap is attached, with the exception of the extreme aft positions in which some trailing- 
edge separation may be present (and desirable). In this region, a radial basis function [18] 
(described in Appendix D.l) serves as the approximation method, which yields the surrogate 
Si(p). In f> 3 . the flow over the flap is fully separated, and a second radial basis function 
fit serves as the surrogate 5.3 (p). In f h. the resolution of the construction points is not 
sufficient to determine the precise location of the separation line. In this region, a simple 
linear triangulation between 5i(p) and S$( p) is used as the surrogate. 5 2 (p). The error 
function. g{ p). is set to unity in and and <?(p) = 50 in reflecting our uncertainty 
in regard to the location of the separation line and, hence, the lack of confidence in the 
quality of the surrogate in this region of the input space. A three-dimensional surface plot 
of the surrogate is shown in Figure 6-5. 

To validate the lift coefficient surrogate, a set of random input points in is first 
selected and the experiment is conducted at each of these points to form the validation 
sample set X va . The input points are confined to the design space Q described in the 
previous paragraph and shown in Figure 6-4. Because the construction data were obtained 
simultaneously with the validation data, there was no expectation in regard to those regions 
of the input space that would be of most interest; thus, a uniform probability density 
function p(p) was used for the selection of the validation points. The number of points 
budgeted for the validation was N va = 45 and, using the relationship in Equation (6.7). 

= 0.03 and £2 — 0.25 were set. If the form of the surrogate had been known prior to 
taking the validation data, the design space could have been restricted to a more feasible 
region and a importance function p(p) that would have concentrated validation points close 


O Construction points 


fi = fiiufi 2 u r> 3 
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Figure 6-5: Three-dimensional mesh plot of the lift coefficient surrogate S{ p). 

to potential designs could have been used. Both of these changes would have improved the 
predictability results for the designs. The scaled model prediction error computed according 
to Equation (6.8) is U = .0482. The maximum un-scaled error does in fact occur in fh_> as 
was presupposed and has a value of 0.4824. If r) ( p ) = 1 had been used everywhere (instead 
of as described above), the model prediction error would have been approximately one order 
of magnitude larger, and would have overwhelmed the results. 

The surrogate just described and the related validation results serve for all of the designs 
discussed in the remainder of this paper. One primary advantage to using the surrogate 
approach is the fact that no additional experimental data are required to bound the errors 
of future designs that are pursued with the surrogate. This characteristic, combined with 
the negligible computational time required for each surrogate evaluation, yields a highly 
flexible design environment that does not sacrifice predictability. 

6.4.2 Single-Point Design, Surrogate Maximization 

For the first study, single-point design is pursued that maximizes the surrogate output. 
The parameter A is set sufficiently large in Equation (6.2) and the resulting function is 
minimized. To accomplish the optimization, an unconstrained quasi-Newton optimizer that 
is included in the optimization toolbox of Matlab [47] is used. The resulting surrogate- 
based optimizer is located at p* = ( x*,y m ) = (.997. .036). and the surrogate-predicted lift 
coefficient value at this point is <S(p‘) = 3.388. The optimizer was started with an initial 
guess at po = (.987, .033) and required 44 surrogate evaluations to arrive at p*. Because 
the surrogate is inexpensive to evaluate (and because there are only two inputs and the 
results can be easily visualized graphically), it can be verified that a surrogate-predicted, 
global maximum is indeed achieved. This verification is more difficult in a purely on-line 
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Figure 6-6: The surrogate-predicted optimizer, p*. and the associated prediction neighbor- 
hood, P(p*,£i). 

optimization setting without restarting the optimizer at multiple initial points po until there 
is sufficient confidence that a global maximum has been obtained. 

Finally, the sensitivity minimizing metric A(a. b) = |<5(a) — *S(b)| in Equation (6.20) 
is selected and the a posteriori error analysis for a single-point design is performed. The 
prediction neighborhood P(p*,£i) is constructed around p* and the surrogate sensitivity 
parameter S — .0328 is found. The optimal point p* and the associated prediction neigh- 
borhood ^(p'.ei) are plotted in Figure 6-6. The resulting predictability statement reads as 
follows: with confidence level greater than .75, a region T C ^(p'bei) of nonzero measure 
exists such that for all p' 6 T 

|S(p')-5(p*)|<e(p*), (6.21) 

where 

e(p*) < Ug{p*.£i) + 6 = .0810 . (6.22) 

The predictability is relatively good with respect to the surrogate-predicted maximum lift 
coefficient, but quite poor with respect to the range of lift coefficients of interest (i.e., 
corresponding to flap positions in fij). 

6.4.3 Multiple- Target Designs 

For the second design study, a multiple-target design is pursued. The motivation for such 
a study might be an interest in examining the lift coefficient at more than one point of 
the deployment of the flap. Specifically, the goal is to to obtain two target lift coefficients: 
A 1 = 3.31 and A 2 = 3.25. Isocontours of the surrogate indicate that a locus of points in fi 
exists for each target that exactly satisfies the design goals. One input point for each design 
is arbitrarily selected: p 1 = (xl 1 !, yt 1 !) = (.987, .033) and p 2 = (x[ 2 ],yl 2 l) = (.979, .033). 
Around each optimizer, a prediction neighborhood chosen from the family of ellipses that 
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Figure 6-7: The surrogate-predicted optimizers and the associated prediction neighbor- 
hoods (shaded). 

have area equal to £\ and centered at p 9 is constructed. The neighborhoods are oriented 
such that they minimize surrogate sensitivity d(p 9 .ei). The optimizers and associated 
prediction neighborhoods are plotted in Figure 6-7. 

For each of the designs ( q = 1.2). it can be stated with confidence level greater than 
0.75 that a region T 9 C V{p q .Si) of nonzero measure exists such that for all p' G T 9 

|5(p 9 ) -S(p')| <e(p 9 ). (6.23) 

where 

efp 1 ) = U + <5(p‘ . £i ) = .0482 + .0198 = .0680. (6.24) 


and 

e(p 2 ) = U + <>(p' 2 .£i) = .0482 + .0201 = .0683. (6.25) 

The above bounds jointly hold on each design. A slightly sharper bound on the average 
error of the two designs is obtained: 

^(p 1 ) + e(p 2 )] < U + .0149 = .0631. (6.26) 

The increased sharpness results from an analysis of the worst-case distribution of the unchar- 
acterized region between the two prediction neighborhoods. Because of the low sensitivity 
of the surrogate in each of the prediction neighborhoods relative to model prediction error 
U, the improvement is slight. 

6.4.4 Comparison with Direct Insertion 

Cases at identical flow conditions have not been examined with both on-line (the method 
of steepest ascent) and off-line (the surrogate approach) optimization methods. However, 
rough comparisons of the resource requirements are of presented and are of guarded use. 
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The on-line results have been reported in an earlier paper by Landman and Britcher [45]. 
In that effort, they found the optimizer to be very robust (successful in 6 out of 6 attempts) 
and insensitive to the initial guess. For each case, they started the optimizer at in initial flap 
position with a low C\ value and obtained a final value within approximately 0.7 percent 
of the maximum C[ value in approximately 20 optimizer steps, requiring approximately 60 
experimental data points (3 points per step). With the surrogate method, 119 points were 
required to construct the surrogate and an additional 45 were used for the validation, for a 
total of 164 experimental data points. For the maximization problem, the a posteriori error 
bound was 2.4 percent of the maximum surrogate value. 

While the surrogate approach seems to compare unfavorably to the on-line method, 
several subtleties lie in its favor. First, for designs chosen with the validated surrogate in 
the future (e.g., the multiple-target design examined in this paper), similar error bounds still 
apply and do not require additional experimental data. In contrast, the on-line approach 
would require additional experimental results. Second, a total of 60 evaluations to obtain 
an optimal point with the on-line method can be deceptive: to be assured that the result is 
indeed optimal, additional information is required. The additional information for the study 
cited was in the form of contour plots of a matrix of data. If visualization is not possible, 
a number of optimizer restarts would be required to be assured of a global optimal. Third, 
in cases for which the objective function is less forgiving, restarts of the on-line optimizer 
would be unavoidable, which would further increase the required experimental data to a level 
surpassing that of the surrogate approach. Finally, the obvious difficulty in pursuing on-line 
optimization is related to the ultimate application: if the intent is to incorporate the data 
as a portion of a larger optimization study, no alternative is available other than to store 
the experimental data for later use and extract with some form of an approximation. If one 
is restricted to a purely experimental setting, then the ability to quickly, and automatically, 
find optimal operating points with the on-line optimizer is highly advantageous. 
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Chapter 7 


Concluding Remarks 


The work presented in this thesis represents several significant contributions. First, the 
level -set based geometry description provides a unified means to describe a family of 
shapes that encompasses multiple topologies. Second, se%-eral theoretical extensions to the 
nonparametric-validated surrogate framework have been derived and demonstrated. Third, 
the synergy that results from the marriage between the surrogate framework and a Pareto 
formulation of the multicriteria design problem have been documented. 

The level-set based geometry definition presented in Section 2.5 represents a way to 
describe a family of shapes that encompasses multiple topologies and that is applicable 
to surrogate-based optimization approaches. The shape description is consistent regard- 
less of the topology of the geometry (single or two bodies) and requires no additional 
logic for the different topologies. The critical characteristic of the method for its appli- 
cation to surrogate-based optimization is that, as shown empirically in Section 2.5.4, the 
input-output functions are continuous across the topology change. Continuity is required 
if interpolatory surrogate models are to be readily fitted through a collection of input- 
output points. The family of shapes defined by the level-set based geometry description 
includes cylindrical eddy-promoter configurations that have been examined previously, al- 
lowing improvements gained by a richer family of shapes to be identified. Additionally, 
the method uses a function superposition technique to define the geometries. It therefore 
can be extended in a straightforward manner by the superposition of additional functions 
to define richer families of geometries with more complex single body configurations, and 
configurations with more than two distinct bodies. 

The surrogate framework has been extended in a number of ways. The formulation 
given in the thesis is applicable to a general, two performance metric, two output design 
problem. The performance metrics can be explicit functions of some or all of the inputs as 
well as of both outputs. The input variables that enter into the performance metrics only 
explicitly (i.e. the inputs that do not enter into the outputs) are handled independently from 
the modeled inputs (the inputs that enter through the outputs) so that as much analytic 
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information as possible is used. By validating the outputs instead of the performance 
metrics, flexibility is afforded the designer in that the validated outputs can be used any 
time in the future for other design problems using different performance metrics. If the 
performance metrics are validated, the designer is restricted to those functions for any 
future design studies. 

A new formulation for the a posteriori error analysis, proximal candidates, has been 
derived. Previously, the proximal region error bounds similar to those given in Sections 
4.4.1 and 5.4.1 were used. The proximal region error analysis bounds the truth values 
on points in a finite-sized subset of the prediction neighborhood. The proximal candidate 
error analysis bounds presented in Sections 4.4.2 and 5.4.2 apply to a randomly drawn input 
points selected from the prediction neighborhood according to the scaled validation sampling 
density. While the proximal region analysis provides a sense of stability to the results (not 
just one. but many designs near the surrogate predicted optimizer satisfy the bounds), the 
proximal candidate analysis gives a probabilistic bound on the truth performance for an 
actual design. If the designer determines the proximal candidate bounds to be satisfactory, 
the design process can be considered a success. In general, the proximal candidate bounds 
are likely to be sharper than those for the proximal region analysis. 

Finally, by combining the surrogate optimization framework with a Pareto analysis for 
a two-criteria design problem, several benefits are realized. First, the Pareto analysis pro- 
vides a great deal of information to the designer as to trade-offs and the results that is gives 
can be rapidly assessed. Little generality, in terms of applicability to other multicriteria 
formulations, is lost as was demonstrated in Section 2.4.2. Second, a full Pareto analysis 
is only possible under the best of circumstances. To be practical, the cost each evaluation 
of the performance metrics must be trivial. This qualification almost never applies for the 
truth (the simulation or the experiment) necessitating a surrogate approach. The combined 
surrogate-Pareto framework takes advantage of the inexpensive surrogates by using them to 
screen out uninteresting (non-Pareto-optimal) regions of the design space. The remaining 
space of surrogate-predicted. Pareto-optimal designs is used for any future design studies. 
Third, and finally, the Pareto analysis improves that predictability of the surrogate frame- 
work, error analysis. The scalarization procedure reduced the effective dimension of the 
design space from the number of inputs ( M ), to one less than the number of outputs. For 
the two-criteria eddy-promoter problem, the surrogates are validated for the predictability 
analysis over only a single dimension. This greatly increases the predictability of the results 
that suffers rapidly with increasing input dimension as demonstrated in Section 4.5. 

The primary weakness of the surrogate framework is that the sharpness of the validation 
results suffers greatly as the input dimensionality increases. This weakness has been de- 
scribed in Section 4.5. For problems in which there are several performance metrics and in 
which there is a large number of inputs, it may be more practical to formulate the problem in 


136 



a Pareto sense. The Pareto formulation reduces the effective dimension of the problem and. 
provided there are only a few performance metrics, the predictability is greatly improved as 
was shown in Chapter 5. The baseline surrogate framework has been shown in prior work 
[77, 79] to be very efficient for problems with several inputs and a large number of outputs. 
The surrogate-Pareto approach provides a strategy for approaching problems with a large 
number of inputs and only a few outputs. Between the two surrogate frameworks, the only 
problems that are completely out of the question are those with a large number of inputs 
and a large number of outputs. 
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Appendix A 


Plane Channel Heat Exchanger 


In this appendix, the solution to the plane Poiseuille problem is presented. The channel 
problem can be solved exactly, and serves as a comparison to the eddy-promoter exchanger 
results presented in Chapter 5. The plane-channel can also be interpreted in a Pareto sense 
as a (admittedly trivial) trade-off problem with a single design variable t)l = jj t. The 
geometry for the plane channel is shown in Figure A-l. The lower wall of the channel is at 
y — — 1 and the top wall at y = 1. 

The nondimensional flow rate in the plane channel 



V L 

Figure A-l: Plane-Channel Heat Exchanger Geometry. 
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V y) 4 l 2 12 2 


(A.4) 


where the superscript ‘ PP 1 indicates the plane Poiseuille solution. The plane Poiseuille 
solution is steady and also invariant in the x-direction. The pressure forcing term on the 
right hand side of the r-momentum equation required to achieve the flowrate Q = | is 


/ 


PP _ 


Re 


(A.5) 


With the exact solutions given above, the plane Poiseuille outputs, analogous to #o(Pm) 
and V'o(Pm) for the eddy-promoter, can be obtained analytically. The temperature output 
can be written as 

(A.6) 


nP P rp rp 

^0 — 1 w ~ 1 m* 


where T w is the wall temperature and T m is the mixed mean temperature 

1 


Tm = Q 1 1 uPP ( y W PP ( y ) dy - 


(A. 7) 


The over-bar in A.6 means to average spatially over the periodicity cell length l = 6.666. 
The channel flow result for the temperature output is 


aPP _ TfT _ 7fT- _ ^3 _ _39_ _ 26 
lm — 16 560 ~ 35 

The pressure output ipQ P is identically given bv Equation A.5 


(A.8) 




Re 


(A.9) 


The performance metrics can be formed by substituting the exact outputs in Equations 
A.8 and A.9 into the performance metric expressions in Equation 2.97 and 2.98. The 
resulting, exact performance metrics for the plane-channel heat exchanger are 


Q PP (vl) = log 10 
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(A.10) 


* PP (riL) 


lo §io [l>o P Re 3 vl 

logio 


) [ 2 Re 2 r,l 


(A. 11) 


For the problems studied here, the Reynolds number has been fixed to Re ~ 250. The 
only remaining independent variable for the heat exchanger is the inverse height parameter 
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T) L . In this case, the Pareto-optimal output manifold dA pp is identically the achievable 
set A pp . It is obtained by evaluating equations A. 10 and A. 11 for a range of r\i € y = 
[VLmin^Lmax]- This has been done for rj Lmin = 20.0, VLmax = 500.0 and the results are 
plotted in Figure A-2. At the top of Figure A-2. the Pareto output curve with the pumping 
power metric 'F PP {t]l ) plotted versus the temperature performance metric 0 PP is given. 
Directly below the Pareto curve is the inverse height parameter tjl plotted versus 0 PP . 

In Chapter5, the plane-channel heat exchanger is used as a reference case for evalu- 
ation of the eddy-promoter heat exchanger performace. Included on the Pareto-output 
plots presented for the eddy-promoter heat exhanger is the plane-channel Pareto curve. 
This curve has been obtained by by solving Equation A. 10 for t)l over the range of tem- 
perature peformance metric values achieved by the eddy-promoter heat echanger. This r//, 
result is the inverse height value required by the plane-channel exchanger to achieve the 
given temperature performance which is then substituted into Equation A. 11 to obtain the 
plane-channel exchanger pressure metric value that corresponds to the same temperature 
performance as the eddy-promoter heat exchanger. 
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Appendix B 


Surrogate Proofs 


B.l Baseline Surrogate Validation Proof 


The relationship between the validation sample size theorem given in Equation 4.5 and the 
resulting validation statement given in Equation 4.8 can be proved using order statistics 
[15]. The proof of this result has been given previously in [36. 37. 77] and is included here 
for convenience. The proof is given for the most simple, single output function case, without 
error scaling functions is given in detail first. Finally, the generalization to the two output, 
scaled error analysis cases used throughout the thesis is given at the end of this appendix. 

It is fist assumed that there is an input-output function <S(p) : 51 — > 1R and a cor- 
responding surrogate S{ p) : Q -> M to <S(p). No assumptions are made as to the form, 
smoothness, or continuity of either «S(p) or «S(p) and. although in practice it is preferred 
that 5(p) approximate 5(p) as well as possible, the results are valid regardless of the quality 
of approximation. The model prediction error function £(p) : f> -» IR is defined as 

£(p) = \S(p)-S(p)\. (B.l) 

A function Z(x) : [0. oc) -¥ [0. 1] is introduced and is defined as 

Z(x)=n p ({pen\£(p)>x}). (B.2) 

The function Z(x) gives the measure of the subset of the input domain 9 for which £(p) is 
strictly greater than x G [0. oc). For situations in which there are regions of finite measure. 
9. c C 0. for which £(p) is constant, {£(p) = x c . Vp € Q c }. it is necessary to define the jump 
in Z(x) as 

lim(Z(a- c — y) — 2{x c )) — fi p (9. c ). (B.3) 

y\0 
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Figure B-l: Plot of the measure function Z(x) versus the surrogate error level x . 


The model prediction error U is computed as 

U= max (B.4) 

where the input points P z . . . , P\ are drawn randomly according to the probability density 
function p(p) : Q — t JR. P; ^ p(p). The random, uncharacterized region is defined as 

T = {p G | £(p) >U). (B.5) 


Next, a random variable Z G [0, 1] which is the measure of the random set T is introduced 
and defined as 

Z = M T), (B.6) 

and the goal is to determine the cumulative distribution function for Z. fz{z) = Pr{Z < z}, 
from which the desired validation statement will follow directly. 

If a variable x z is introduced for z € [0. 1] and defined as 


x- = min x. 

{x<=[0 r 3C)|Z(x)<--} 

it follows that 

Pr{Z(U) < z} = Pr{U > x : }. 
Introducing a set V C Q defined as 


V = {p e | £(p) > Xr}. 


(B-7) 


(B.8) 


(B.9) 
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then the measure of the set V, is given as 


H P (V) = inf Z{x) (B.10) 

{*<*=} 

and can be bounded as 

^p(P) > ( B -ll) 

From the definition for the model prediction error, the following statement can be made 

Pr{U > x z } = Pr{3j e (l,...,iV) | Pj € V} 

= l-Pr{Pj en\V,Vj € .N)}. (B.12) 

By making use of the relationship that the validation points P, are i.i.d.. the probability 
that they will all lie outside of V can be evaluated as 

Pr{Pj € Q\V.yj <= = (Pr{P 1 6f!\P}) V 

= ( M „(0\P)V v 

- (l- Mp (P)) v . (B.13) 

From B.10. it is known that (1 - p P (V )) < (1 — z). which gives 

Pr{Pj 6 n\V. V; 6 (1 .V)} < (1 - =) X ■ (B.14) 

and the cumulative distribution function F : {z) can then be bounded by FMz) as 

F-.{z)>F : {z). (B.15) 

where 

F.{z) = 1 - (1 -*)- v (B.16) 

For cases in which Z(x) is continuous. F z {z) = F : {z). Note that, even with the strict 
inequality. Pr{Z < z} satisfies 

Pr{Z < z) > Pr{F z {z)} (B.17) 

although Pr{Z < z) > Pr{F z {z)} can be false for Z(x) discontinous. 

Recognizing that ei can be substitued for z in Equation B.16 and that £2 = 1 
and the solving for the sample size N, the validation sampling theorem follows 


ln(l - £i) r 


-F : (z). 

(B.18) 
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where \z ] is the smallest integer that is greater than z. 

The extension to the two-output analysis cases examined in the body of this thesis 
is straighforward. By redefining the error function £(p) to be identical to the function 
validated in Equation 4.6 


£ (p) = max 


|^(Pm,)-^l(Pm,)| |#(Pmi)-fc(Pmi) 


(Pmi) ffdij(Prai) / 

the remainder of the analysis follows exactly as has been outlined above. 


(B.19) 


B.2 Predictability Analysis: Proximal Region 

The proof of the predictability results given in Section 4.4 is based on the validation result 
in Equation 4.8. For a point p' m to exist such that at least one of the outputs bounds 

Id <<Pl(p'm) <«0P (B.20) 

^o> < </>f(p',J < u 02- (B.21) 

are not valid, the point would have to be in T. This follows from the result that if p'^ E 
n' m \Y, then 4.8 holds, and the bounds given in Equations 4.16 — 4.19 hold as well. 

The prediction region V* : is selected by the designer. The worst case scenario in terms 
of the error bounds is for the prediction to fully contain the uncharacterized region T. If 
this is the case, then the region T f E V* x will exist given the strict inequality on the size 
of T in Equation 4.8. This result highlights the necessity of having a prediction region V* 
that is of at least p- measure S\. 

The proof of the performance metric bounds in Equations 4.22 — 4.27 is similar. Because 
the output have been shown to be correctly bounded by 4.16 — 4.21 for p' m E then for 
points p' E P, the metric bounds will hold as well. 


B.3 Predictability Analysis: Proximal Candidates 

The proof of the proximal candidate predictability result begins by considering the the error 
function £(p) : ft — > 1R defined as 


£(P) = |S(p)-5(p)|. (B.22) 

The validation sample set is formed by selecting N points. Pi. . . . Py. randomly according 
to probability density function p(p). The validation sample set X va is 

= {(P 1 ,5(P 1 )), . . . , (P jV!< S(P lV )). }, P ; ~ p( p). (B.23) 
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and the model prediction error is computed as 


U = max £(P,)= max |5(P,) - 5(P ; )|. (B.24) 

ie{l A'} 

Next, the random variable z is introduced and is 

z = J_p(p)dp (B.25) 

where 

E = {p e n | £(p) < U). (B.26) 

The cumulative distribution function for z. Fj(z). is 

F-(z) = z s . (B.27) 

which can be easily shown by order statistics [15]. The complement to z, z = 1 — z is 
introduced and it follows that it has cumulative distribution function 

F : (z) = 1 - Fs[z) = 1 - (1 - ^)‘ v . (B.28) 


The random variable z is the p-measure of the uncharacterized region T. Next, consider 
a prediction region V a , of p - measure a. The expected fraction of designs V a for which 
£(p) < U can be bounded by 

a> Jo f:{z) dz (B,29) 

where f : {z) is the probability density function for z and is 

f : (z) = = iV(l - z) x ~ l . (B.30) 


The integration in Equation B.29 gives the conditional expectation for the uncharacterized 
region completely inside of the prediction neighborhood. Substituting B.30 into B.29 gives 


a > 


f (!-?)",* -O'' 1 * 


1 

a(iV + 1) 



(B .3 1 ) 
(B.32) 


Finally, recognizing that a = 1 - e c , the expression in Equation 4.37 and given in Section 
4.4.2 is obtained 


£ c 


1 

a(N + 1) 


(i-d-cr)^ 1 ). 


(B.33) 
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P2 A 



Figure B-2: Schematic of the input domain fh prediction region 7V and uncharacterized 
region T which is contained completely inside of TV 


B.4 Optimality Analysis 


For a sequence of identically distributed, randomly selected inputs Pi,....P ; , is drawn 
according to the probability density function p(p m ). The random variable 

L = min j such that $i(P j) < $i m j n , and ^2(Pj) < $2min! (B.34) 

is introduced. It can then be shown that [36], for the validation sample size iV, that 

Pr{L > 1} > (B-35) 

It follows that, for ei €]0, 1[, that if m is set such that 

m < N (J- - l) , (B.36) 

that the probability of drawing a sequence of random input vectors Pi,... . P m according 
to p{p m ) and finding a point such that both 


$ 1 (0 1 (P mj ),^ 2 (P mj ),P J ) < $ lmin! (B.37) 

< ^2min: (B-38) 
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is less than l - El- The lower bounds on the performance metrics are given by 


^lmin 

min 

$(zi,z 2 .p'), 

(B.39) 


{p'€n',2!€2° pt , 



< ^ > 2min 

min 

$(zi. 2 - 2 -p')- 

(B.40) 


where Z\ and Z-i are 

ZT = {z | vp' m € n' m . [Mp'J - Ug 0l ( p'J] < z < ifdp'J + Ug 0l (p'm)]}- (B.41) 

Z° pl = { 2 I V P ;„ € fi' m , [Mp m) - Ug 02 {p' m )) < z < [fc(p' m ) + Ug 0l (p' m )]}. (B.4‘2) 

For a point p' to have truth peformance metrics such that the lower bounds and 

$2min in Equations B.37 and B.38 do not hold, the point must lie in T. If p' € T, then the 

output bounds 


[<Al(Pm) - U 9oi( Pm)] < ^l(Pm) < [<MPm) + ^Oi(Pm)] (B.43) 

[ f p2 (P;n ) ~ U<)o>( Pm)] < <h(Pm) < [<MPm) + C<7o ; (Pm)] < B ' 44 ) 


embedded in Equations B.41 and B.42 will not hold, and the performance metrics evaluated 
at p' can satisfy Equations B.37 and B.38. Note that if the point p' is in T, there is no 
guarantee that the performance metrics will satisfy B.37 and B.3S. only that the possibility 
exists, and the confidence level 1 - Zi represents an upper bound. 


A random variable z = p p {Y) is defined as the p - measure of the uncharacterized region 

T. It follows that i 

Pr{L > 1} = f Pr{L>l\z}dF z ■ (B.45) 

Jo 

It is assumed that the point p ; E TZ C T and that the measure of the subset 7 Z is given as 
h(z) — f.ip( 7Z) < z. The probability that of the l randomly drawn points that none of them 
lie in 7Z is given by 

Pr{P* <LH{z). J = l.....l} = (l-h(z)) 1 . (B.46) 

which is then substituted into Equation B.45 to obtain 


Pr{L>l } = [ l Pr{L > l \ z}dF z (B.47) 

Jo 

= [ (1 -h(z))‘dF z (B.48) 

Jo 

> f (1- z) l dF z (B.49) 

Jo 

> [\l-z) l dF Z - (B.50) 

Jo 
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The step in B.49 follows from the relationship that h(z) < z . and B.50 uses the results given 
in Equations B.15 and B.16* Finally, substituting Equation B.16 into Equation B.50 and 
evaluating the integral gives 


Pr{L>l } > f (1 -z)‘N{l -z)‘ w -\ 

Jo 

N 

n + t 


(B.51) 

(B.52) 


At this point, all that is left is to recognize that Pr{L > 1} = e L and to substitute into 
Equation B.52 to get the expression in Equation B.36. 
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Appendix C 


Numerical Methods 


C.l Gaussian Quadrature 

The general formula for the Gaussian quadrature used to evaluate the elemental matrices is 
detailed below. To approximate the integral of a function r/(x) : Q kc — > 1R over a triangular 
element. f>* e . the following quadrature rule is used 

ff(x) dx = ( C1 ) 

1=1 

where .4°*' is the area of the element, sp is the quadrature weight at that point, and x(£,) 
is the quadrature point position in terms of the barycentric coordinates £. A Gaussian 
quadrature integration scheme that can exactly integrate a 5'^-order polynomial requires 
M q = 7 quadrature points. The weights ct, and weight -points £, = (£iof'2,- 6^ * = 
1. . . . .Af 9 are listed in Table C.l [39]. 



i 


fli 

bi 

$3i 

i 

0.225 

1/3 

1/3 

1/3 

2 

0.1323941527 

a 

b 

b 

3 

0.1323941527 

b 

a 

b 

4 

0.1323941527 

b 

b 

a 

5 

0.1259391805 

c 

d 

d 

6 

0.1259391805 

d 

c 

d 

7 

0.1259391805 

d 

d 

c 


Table C.l: Fifth-order. Gaussian quadrature points, M q = 7: a = 0.0597158717, b — 
0.4701420641, c = 0.7974269853, d = 0.1012865073. 
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n — 2 n — 1 n n + 1 

Figure C-l: Uniform-in-time difference mesh 


C.2 Variable Time-Step, 3 rd — Order Adams-Bashforth 

C.2.1 Fixed Time-Step Schemes 

The Adams-Bashforth difference schemes are a family of explicit, multipoint integration 
methods [11, 27]. In general, a fc th - order Adams-Bashforth scheme requires k function 
evaluations at time-steps n.n — 1, . . . , n — k 4- 1. The general problem 

u(t) = /(u(f). t) (C.2) 

will be central to the discussion of this section. 

The uniform, finite difference grid upon which the derivation is based is given in Figure 
C-l. To derive the 3 rd -order Adams-Bashforth scheme for fixed time-step. A t. the general 


form of the resulting scheme is assumed to be 

Wnf 1 = Un + At(Pofn + Plfn- 1 + Pifn)- (C.3) 

where the u rn = u(mAt), f m = f(u(mAt), mAf). and /3 t . i = 0.....2 are the unknown 


coefficients that must be determined to yield a 3 rd -order scheme. Next, each term is ex- 
panded about time-step n with a Taylor series, noting that from Equation C.2 u = /', and 
v! = /". The Taylor series are substituted into Equation C.3, and the expressions for the 
coefficients are solved that eliminate 2 nd -order and lower error terms. The result is the 
familiar 3 rd -order, fixed A2, Adams-Bashforth scheme [27]. 

A t 

^n+i ~ Un d - qTr(23/ n — 16fn—i -h 5/ n ). (C.4) 


Because, the 3 rd -order Adams-Bashforth scheme is fully explicit and has a very good 
combination of accuracy and stability properties, it is an attractive choice to advance the 
nonlinear convective part of the Navier-Stokes equations. The intersection of the region of 
stability with the Imaginary axis is 0.723. For the pure convection equation u t = uu x . the 
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Figure C-2: Variable time step size difference mesh 


stability limit is thus given as 

A t < 0.723 

max (u) 


The stability limit is known as the Courant condition. 


(C.5) 


C.2.2 Variable Time-Step Schemes 


For problems in which the convective velocity varies significantly over time, a fixed time- 
step scheme is impractical. A solver constrained by a fixed step would require that the user 
be able to anticipate maximum convective velocity a priori so that the time-step could be 
set at the beginning of the run according to C.5. Not only is it difficult to anticipate the 
maximum convective velocity, for problems in which it varies significantly, such a restriction 
greatly increases the required number of steps and hence, the required computational time 
to obtain a solution. For these reasons, a variable time-step size version of C.4 is derived. 
In the following section, an Adams-Bashforth scheme consistent with the 2 nd -order time 
advancement is derived. 

To proceed, the non-uniform finite difference grid shown in Figure C-2 is considered. 
Using the method of undetermined coefficients, the terms in Equation C.3 are first expanded 
about time-step n with a Taylor series; 


*n-r 1 


n 2 o 3 

. f , u ft | u „ ttl | r\f n ^ \ 
= u n + au n + —u n + Y n + ‘ 


(C.6) 


-l = u n 


, , ^ b 2 „ 

bu n + ~ 2 U n 


m 

6 U " 


0(b *) 


(C.7) 


In = < (C.8) 

Sn- 1 = < - K + + 0(b<) (C.9) 

in-2 = < - : + jK - p;: + o(c 4 ) (c.io> 

where, in the interest of brevity, a — A t n . b = Af„_i. and c = Af„_i + A t n -2- If Equa- 
tions C.6— C. 10 are substituted into C.3, and like terms are collected, the following set of 
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equations that must be satisfied to obtain a 3 rd -order accurate scheme are obtained; 


’ll 1 


’ /3o ' 


a 

0 —b — c 


Ai 


a 2 

2 

b 2 c 2 




a 3 

y y . 


. & . 


- If - 


(C.ll) 


The matrix Equation C.ll is solved for /?o, P\, and /? 2 . Substituting for a. b. and c, and 
simplifying gives the variable time-step size, 3 rd -order Adams-Bashforth scheme [28] 


2 

u n-r 1 ~ 4 ^ " Pi fn-i 

i=0 


(C. 12) 


where. 


Po ~ 
Pi = 

Pi = 


Af„ 

|" 12A^_i(At ri _i 4* A* n _ 2 ) + 6At ri (2At n _i 4- A£„_ 2 ) 4- 4A^ 

12 

At n 

[ At n _i(At„_ 

[6At n (A< n _i 4- At n _ 2 ) 4- 4Af 2 ] 

1 + A£„_ 2 ) 

^ <1 

[ Af n _iAt n _ 2 

6At ri At„_i + 4 At 2 

. 

12 

A^ n _ 2 (At rt ,i 4- A f /j _ 2 ) 



(C.13) 
(C.14) 
(C. 15) 


As expected, when — Af /Z _i = A t n = A^,- the scltem^-wi Equation C.4 is recovered. 

The derivation of the 3 rd -order Adams-Bashforth scheme consistent with the 2 nd -order. 
backward difference, time advancement of the velocity is similar to the derivation just 
described. In this case, the resulting scheme must satisfy the difference form 


ot\u n + 1 + a 2 u n - a\u n -\ = Pof n 4- Pif n - 1 4- /?2/n-2* (C.16) 


where 


ai 


2A t n 4- A^_i A t n 4- A£ n _i A t n 

■ ■ — ■ - , Q; 9 — 1 ~ 

A£n(A£ n + A i n _j) A£ n A£ n _i A£ n _i(A£ n + At n _i) 


(C.17) 


Again, the Taylor series substitution into Equation C.16 is made and solved for /3;, 
i = 0. . . . , 2 that eliminate all errors that are 2 nd -order and lower. The resulting coefficients 
are 


A) 


Ai 


A tn 
12 

A t n 
~12 


+ Atn- 2 ) + 6At fl (2At„_i + A£ n _2) 4- 4At 2 
At n _i(At„_i + A<„_2) 

6At ft (At n -i + At n - 2 ) + 4A^ ' 

At„_iAt n _9 


(C.18) 
(C. 19) 
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02 = 


At n | 6At n At n -i + 4At 2 n 


12 Ai n _ 2 (At n _i + At n - 2 ) 


(C.20) 


For the fixed time-step case in which At n -2 — At n -\ = A t n = A£, the coefficients reduce 


_ A^ n _ Af n _ 

A) - 12 - .1 12 ’ 12 ' 


(C.21) 
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Appendix D 


Surrogate Models 

D.l Radial Basis Function 

We use a radial basis function model as a surrogate. Radial basis functions are a well 
known class of frictions used in a number of applications [18]. We begin the description of 
the surrogate used by assuming that we want to construct a surrogate model valid in the 
design space 0. G 2R U for the general function /( p) : Q -> 2R. /( p) ~ / (p)- The radial 
basis function surrogate has the general form 

yco 

/( p) = a i&( r i(p)) + ^«(p) (D-l) 

1=1 

where (^(r^p)) is the radial basis function and P m (p) is the polynomial of degree m-1. For 
a d-dimension input space, we use the notation p = (pi . . . . - Pd } and p™ = {p™ • • • ■ • P c t 0 d)- 
Given this notation, the radius r,(p) is simply the Euclidean distance from construction 
point i 

n(p) = L(Pfj-Pj) 2 

J = l 

For the surrogate used in the work reported here, we have used the radial basis function 
cp(r) = r 2 ( m-1 Mogr and have set m = 2. The coefficients a l of the general form (D.l) are 
determined from 

,\ co d 

y <w,(r,( P n) + E w p") = /(pf ). * = i jv" (° si 

1=1 j= 1 

x co 

pf°) — J = ( D - 4 ) 

1=1 

The polynomial contribution for m = 2 is q(p™) = {l.p^. . . . .p^} T - 
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Appeals to the simulation provide /( pj°). The solution for the unknown coefficients 
a = {gi , . . . . a^v co } 7 ' and b = {b\< . . . , 6^+1}^ is the solution to the following matrix equation 


A rbf E rbf a _ f 

^rbf 0 b 0 


(D.5) 


where f = {/( pf 3 ), . . . , (p*y co )} T is the solution vector at the construction points. There is 
a very real possibility that the matrix equation (D.5) could be ill-conditioned [18], but for 
all of the models generated here for a modest number of construction points (on the order of 
a couple of hundred), we have not experienced any difficulties with the required inversion. 


D.2 Orthogonal Arrays 

A critical consideration in approximate optimization approaches is how to best select con- 
struction points from the input domain Q in order to extract as much information as possible 
without having to sample at an inordinate number of points. The method used to select 
points for the output surrogates constructed in Sections 4.1 and 5.1.1 is based on the work 
of Bose [8. 9]. The methods have formulated and coded into C-routines by Owen [58], and 
are part of a very useful package for generating experimental designs. 

The array that has been generated for the eddy-promoter surrogate construction point 
selection has k — 3 columns, strength t = 2. and q = 11 levels [58]. An orthogonal array of 
strength t is a matrix with n rows and k columns such that for any ti x t sub-matrix, the q l 
possible rows occur the same number of times. The strength k = 2 array with q = 11 levels 
has N orth = q k = 121 rows (input points). The rows were permuted and the construction 
points were selected from the N ortk = 121, randomly permuted, orthogonal array points. 
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